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PREFACE 


This Memorandum Report consists of a compilation of three individual 
reports, of increasing complexity, describing investigations of formation flight 
of spacecraft in the vicinity of the L-i Sun-Earth libration point. The indi- 
vidual reports form the following parts of this compilation: 

• Introduction to the relative motion of spacecraft about the Sun-Earth 
L 2 Point 

• Linear and quadratic modelling and solution of the relative motion 

• Modelling the Perturbations — Elliptical Earth Orbit, Lunar Gravity, 
Solar Radiation Pressure, Thrusters 

The three parts are self-contained, with somewhat, varying notation and 
terminology. 

This work was funded by the Distributed Spacecraft Technology Program 
at NASA's Goddard Space Flight Center. The sponsor was Dr.. Jesse Leitner 
of the Mission Engineering and Systems Analysis Division. 

After fairly significant literature searches, this new work (of Parts 2 and 3) 
is deemed to be unique because it describes the primary perturbations to the 
description of relative motion between nearby spacecraft. The effect of the 
elliptical motion of the Earth about the Sun was verified to be the dominant 
perturbation to the circular restricted three body problem. Contributions 
due to lunar gravity and solar radiation pressure arc seen to have much 
smaller effect. 
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1 Introduction 


The long-term goals of our sponsor are to develop high-fidelity equations 
of motion representing orbiting spacecraft near the Sun-Earth Lo point and 
equations of relative motion between orbiting spacecraft near the L 2 point. 
These equations will be used to develop orbit control schemes for a constel- 
lation of spacecraft near this Lo point. 

This report provides a starting point to the future development of high- 
fidelity equations by first developing equations of motion for the circular 
restricted three-bodv problem. Additionally, equations are derived for ana- 
lytical expressions for the relative motion between spacecraft orbiting near 
the L -2 point. 

We follow the terminology of Hamilton’s thesis [1] as we present four 
sets of equations of motion: three nonlinear sets and a linearized set. Their 
derivations are in separate appendices. These equations were developed to 
examine relative errors of the linearized and quadratic effect models with 
respect to the baseline set full nonlinear equations. 

The four types of models for the circular restricted three-bodv problem: 

• full nonlinear model of spacecraft motion with respect to the system 
barycenter 

• hill nonlinear model of spacecraft motion with respect to the L 2 point 

• linearized model of spacecraft motion with respect to the Lo point 

• quadratic (nonlinear) model spacecraft motion with respect to the L 2 
point. 

We formed a Taylor series expansion of the nonlinear equations about 
the libra tion point. Our linearized and quadratic, models are respectively the 
first- and second-order Taylor series. 

The analytical expressions of relative motion are based on the solutions 
to the linearized equations. 

2 Problem Definition 

The restricted three-body problem is a simplification to the general three- 
body problem. The general three-body problem is the description of orbital 
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motion of throe bodies, of arbitrary mass, in arbitrary orbits. The general 
three-bodv problem has not been solved in closed form. ‘‘Work has focused 
on simplifying the general problem. One special analytical solution - the 
restricted three-body problem has been known since the time of Euler and 
Lagrange.'’ [2] 

“We define our problem as follows: Two bodies revolve around their cen- 
ter of mass in circular orbits under the influence of their mutual gravitational 
attraction and a third body (attracted by the previous two but not influenc- 
ing their motion) moves in the plane defined by the two revolving bodies. 
The restricted problem of three bodies is to describe the motion of this third 
body T [3] 

The two revolving bodies are called the primaries. The masses 774 and 
/77 2 of these bodies are arbitrary and are considered as point masses. The 
mass of the third body, 7773, is infinitesimal and does not influence the motion 
of masses mi and m 2 . I11 our problem, consider the mass to be a space- 
craft. Also consider the sum of the mass of the Earth and its Moon to be 
mo. Assume the Earth and Moon are in circular orbit about their common 
bary center (system's center of mass). In turn, this baryeenter is in circular 
orbit about the Sun, which is mass ni\. These four masses compose our 
system; however, there are two large bodies. The larger of the large bodies 
is the Sun and the smaller of the large bodies located at the Earth-Moon 
baryeenter (774 > ni 2 ). 

We copy portions of Hamilton’s thesis Section 2.1 to define terminology 
and coordinate systems. 

For this preliminary analysis, the orbits of the Moon around the Earth 
and the Earth around the Sun are assumed to be circular, have a constant 
angular speed, and be in the same orbital plane. (Spacecraft motions are not 
restricted to the plane.) When these assumptions are made, the restricted 
three-body problem becomes the more specific circular restricted three-body 
problem. 

The rotating coordinate frame shown in Figure 1 is used for the devel- 
opment of the equations of motion. The origin is the system baryeenter of 
777 1 and m2. The first basis vector, 74 , points toward the smaller body. The 
third basis vector, o :i , is parallel to the orbit normal with origin at the system 
baryeenter. The second basis vector, is defined by the cross product, of <7,3 
and di (a 2 — a 3 x a j). These vectors form an orthogonal coordinate frame. 

The position vector of the spacecraft in this rotating frame is 
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Figure 1: The Three-Body Problem’s Rotating Coordinate System. 


v — A Q-i + Y d 2 H - Zd 3. 

Writing vectors from the large bodies to spacecraft in the rotating frame 
gives 

i'\ — (A H- Di)&i -(- V do + Zd 3 

and 

v 2 — ( X — D 2 )d 1 + Yci2 — Za,i, 

where D\ is the distance from the barycenter to in 1 and D 2 is the distance 
from the barycenter to mo. Note that 


D — D\ + D 2 . 


Use the equation for the definition of center of mass: ni\ D\ = m 2 D 2 to 
calculate the magnitudes of and Z?2 given D: 


D\ 

D 2 


D 


1 + 2 

rr 

D 


1 + 


m\ 
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Define the constant angular speed of the frame about d 6 as 


_ jG(?n 1 + m 2 ) _ jG(niSun + Earth + 

“ V Ip ” V “d* ' 

where G is the gravitational constant. As shown in Appendix A. we 
define the constant angular speed of the frame about a$ as the two-body 
mean motion Cj = ujci$ 

Also define 


}t\ — Grrii — Gmsun 

and 

//*2 (j{)f‘2 — G(lllf? nr f] r A' til \J non) ■ 

With these definitions, we state four sets of equations that describe the 
motion of the spacecraft. The derivation of these equations is shown in 
respective appendices. 

3 Full Nonlinear Equations for Circular Re- 
stricted Three-Body Problem (relative to 
barycenter) 

These equations describe the spacecraft’s motion relative to the system's 
barycenter in a rotating coordinate frame. The frame rotates with a constant 
angular speed. 


X - 2ljY - u 2 X = 
Y + 2ujX - u ?Y = 
Z = 


/U-A + Pi) _ //.,(A - P 2 ) 
|n| a |r a | 3 

/q y_ _ my 
I r 1 1 3 up 

/C| Z fl) Z 


( 1 ) 

( 2 ) 

( 3 ) 


where J/q | and | / 2 1 arc the magnitude's' of the vectors i : \ and r > . respectively. 
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Figure 2: Libration Point Locations. 

The initial conditions to be specified are 

*( 0 ) X(0) 
y(0) Y( o) 

m m. 

These equations are derived in Appendix A. 

By setting all time derivatives in these equations of motion to zero, five 
libration points can be calculated. The location of these points depend on the 
masses and distances of the bodies: however, three points, shown in Figure 2, 
are always collinear with the two large bodies (L i, L 2 , L:$), and the other 
two points form equilateral triangles with the two large bodies (/^ and L h ). 

The motion near the collinear libration points is always unstable due to 
the existence of a positive real root of the characteristic equation for any 
value of fi. An unperturbed object in an orbit around a collinear libration 
point will move away from that point. 

The four initial conditions corresponding to the in-plane motions (X-Y 
plane) cannot be arbitrarily selected. This will be further explained in Sec- 
tion 5. Szebehelv discusses that the selection of initial conditions at the 
collinear libration points show considerable inherent instability. Even when 
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Figure 3: Relative Position with respect to Lo. 


initial conditions are chosen according to the linearized requirements, the 
spacecraft orbits maybe once or twice around L 2 before moving away. Fortu- 
nately, periodic orbits may be obtained by slight modification of the initial 
conditions which follow from the linearized solution. The point is that the 
orbits are very sensitive to changes in the initial conditions. 


4 Full Nonlinear Equations for Circular Re- 
stricted Three-Body Problem (relative to 
a libration point) 

We continue from Hamilton. 

Motion around a collinear libration point is more easily expressed in a 
local coordinate frame, shown in Figure 3. with the origin located at the 
libration point of interest. 

Note the vector addition 

r = r () + p 

where fi, has components (A*o, Vo, Z {) ) to define the location of the libration 
point relative to the system barycenter. The ”0 r subscript refers to any 
one of the libration points. (Vo and Z 0 are both zero for collinear libration 
points.) Vector p has components (A, V', Z) to specify the location of 
relative to the libration point. Using 

A = A'o 4- x (4) 
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Y 

Z 


(5) 

( 6 ) 


= Vo + y 
= Zo + z, 

be reminded that (X, Y\ Z) is the spacecraft position relative to the system 
barycenter. The radial direction, x, is eollinear with dj, going from the Sun 
(or largest body) through the libration point. The cross-track direction, 2 , 
is parallel and in the same direction as S3, the orbit normal. The component 
y is in the direction of d 2 - When circular motion is assumed, the y direction 
is the tangential velocity direction of the libration point around the system 
barycenter. These axes form an orthogonal coordinate system. Motion in 
the x-y plane is referred to as in-plane and any motion in the 2 direction will 
be referred to as out-of-plane. 

The equations below describe the spacecraft’s motion in a local coordinate 
frame relative to the libration point (origin) chosen by selection of ( Xq , Vo, 
Z 0 ). Unique distances corresponding to each of the five libration points exist 
for each set of (Ao, Yq, Z 0 ). 

The equations below are the full, coupled, nonlinear equations of motion 
for the circular restricted three-body problem. The coordinate frame rotates 
with a constant angular speed. 

Equations (1) - (3) yield by simple substitution of Equations (4) - (6): 


x — 2 ujy — u> 2 (X 0 + x) 
y 4* 2lCX — a;- (Vq 4- y) 


Hi(Xq + x 4- D 1 ) /x 2 (Ao -I- x — D 2 ) 

r 3 r 3 

/ x r 2 

lh(Yo 4- y) H 2 ( Eq + y) 

r 3 r 3 

7 1 7 2 

/q(Zo + Z) ^{Zq + z) 


where r\ and r 2 are now defined below as 

x 1 = yj (Ao + x 4- Di) 2 + (Vq + y) 2 + (Zq 4- z) 2 

To = \j (Ao 4- x — D 2) 2 4- (Vo 4- y) 2 + (Zq 4- z) 2 . 

The initial conditions to be specified are 

x(0) i(0) 

3/(0) y(0) 

m m- 


11 


Those equations are derived in Appendix B. 

For the L 2 point, where only A' 0 4 0. the equations further simplify. 

5 Linear Equations for Circular Restricted 
Three-Body Problem 

At a chosen equilibrium point , the nonlinear, full, equations of motion can be 
linearized in order to exploit standard solutions to linear equations. The lin- 
ear equations of motion about the local coordinate frame (origin at eollincar 
libra tion point Ln) are given here: 


r — 2 ll'ij — U xx* — 0 
y + 2ujx — Uyyy — 0 
£ — Uzz z ~ 0 : 


where, evaluated at the Lo collinear libration point, where A' = A' 0 , Y = 0, 
Z = 0 


u yy\l 

Uzz\ l 


'■ + — ^ 

(V + D x f (A'o -D 2 f 

. fM_ M2 

(A'o + D,Y (Xo-DX 

Mi t°2 

(X 0 + DX (A'o -W 


The initial conditions to be specified are 

40) x(0) 

40 ) 2 /( 0 ) 

4 . 0 ) 40 ). 


Tliese equations are derived in Appendix C. 

When evaluated at L 2 Uxx, Uyy-. and Uzz, which were formed from U, 
are constant. Tlie pseudopotential U is the centrifugal plus gravitational 
force potential defined as 


U = 


Y 2 (X 2 + Y 2 ) + ^ + 

2 n 


M 2 

>'2 
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and 


n = y /(X + D 1 )* + Y* + Z* 
r 2 = J{X - Do) 2 + Y 2 + Z 2 . 

The respective double-lettered subscript indicates the second derivative 
of U with respect to that lettered- variable. 


The selection of initial conditions for the first-order linear model 
summarizes the discussions of Szebehelv [3], Farquhar [5], and Wie [4]. The 
in- plane characteristic equation of the linearized equations of motion about 
the collinear equilibrium points is of fourth degree. Solution yields two real 
roots (equal in magnitude, but opposite in sign) and two imaginary roots. 
The in-plane position has a convergent mode (due to the negative real root), 
a divergent mode (due to the positive real root), and an oscillatory mode. 
All three collinear libration points are unstable. 

The out-of-plane characteristic equation of the linearized equations of 
motion about the collinear equilibrium points is of second degree. The two 
eigenvalues are imaginary. Thus, out-of-plane motion is simply oscillatory. 

For the linearized equations of motion, these authors show that a solution 
to the equations can be made to contain only the oscillatory modes with 
proper choice of initial conditions: 

Choose 


m = 


( 7 ) 

m = 

-ku xy x(0) 

( 8 ) 


^,2 _j _ jj 

where k = ^ and oj xv is the nondimensional frequency of the in-plane 
oscillatory mode calculated from 


Mx y 

A 

A 2 


— i + P i + 02 

~ 2 — (U XX + Uyy)/ 2 

— —UxxWyy, 
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which moans that onto initial conditions .v( 0) and ij( 0) arc selected, the 
corresponding initial velocity components cannot be chosen at will. 

Although the choice of the 2 initial conditions is arbitrary. Wie s tates 
that choosing x(0) = 2(0) = 0 and £(0) = —y(0)u z (with u z = \J\Uzz\)*. the 
solution further reduces to a quasi-periodic Lissajous trajectory. We will not 
discuss this trajectory at this time. Wo just list his example initial condition. 

Szebchely continues by noting Equations (7) - (8) do not hold when 
higher-order terms in the function U are retained. Nevertheless. the set of 
initial conditions given by these equations furnish the starting point of an 
iteration which leads to the proper initial condition for the nonlinear case. 
The closer the initial point (;c(0). y( 0)) is to the chosen libration point, 
the better a differential correction scheme will work to furnish the initial 
conditions desired for the establishment of periodic orbits in the nonlinear 
problem. 

The linearization process allows only first-order terms in the coordinates 
and velocities: consequently, meaningful results must be connected with small 
values of these variables. Even when the initial conditions are chosen accord- 
ing to the linearized requirements, the spacecraft perforin only one or two 
orbits before they depart from the area. Again, periodic orbits may be ob- 
tained by slight modification of the initial conditions which follow from the 
linearized solution. The point is that the orbits are very sensitive to changes 
in the initial conditions. 

The linear model is amenable to linear analysis methods, which are more 
tractable than nonlinear solutions. The linear model is most helpful for 
preliminary control analysis as discussed by Hamilton. The nonlinear insta- 
bilities of this system must be restrained to maintain the desired periodic 
orbit motion. 

6 Quadratic Equations for Circular Restricted 
Three-Body Problem 

We went one-step beyond the linear dynamics model of Section 5. By further 
expanding U, we obtained a second-order model for the restricted problem. 
This model should yield trajectories closer to the nonlinear model than does 
the linear model. 


= 2 Ax - \bC2x 2 - y 2 - z 2 ) 

= —Ay + 3 Bxy 
= -Az Wxz 


Ml M2 

fo + ^Di ) 3 (x 0 -D 2 y 

Ml M2 

Co + A ) 4 (A ' 0 - D 2 ) 4 

The initial conditions to be specified are 

t( 0) i(0) 

2/(0) M(0) 

2 (0) i(0). 

These equations are derived in Appendix D. 


x — 2u;y — uAr 
y + — l j 2 y 


where 


A = 


B = 




7 Linear Relative Motion Equations for Cir- 
cular Restricted Three-Body Problem 


Consider the linearized equations, with initial conditions selected such as to 
avoid exciting the divergent mode. The solution then takes the form 


X 


Xq cos f }(t — t 0 ) + ^ sin n(t — to) 

y 

= 

—k. r 0 sin tl(t — to) + yo cos Cl(t — to) 

z . 


Z()cmu(t — to) + sinu;(# — t 0 ) 


where subscript 0 refers to the value at time t = t.Q. This solution can be 
written as a pair of decoupled matrix equations of similar form: 



cos Q(t - to) 

jr, sin Q(t. — to) 

Xq 

—A: sin fi(t — to) 

cos fl(t- — to) 

_ yo , 

COS CO (f - to) 

£ sin co (f - t 0 ) 

^0 

— oj am co(t — to) 

COS io{t — to) 

Zq _ 


( 9 ) 

( 10 ) 
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For notational convenience, lot the vectors x and z refer, respectively, to 
[;r (j\ T and [z i] T . Therefore, Equations (9) and (10) may be written as 


x = A kM (t,L o)x 0 (11) 

z = A WM (t,t Q ) z 0 , (12) 


T coa9.{t -t 0 ) jsinft(f-fo) 
fc ' s! ’’ 11 — AusinQff — hi) cos9(t — t„) 

is the generic state transition matrix. Because of the obvious similarity, the 
relative motion solution may be developed for one of these matrix equations, 
and the result applied to the other. 


7.1 Relative Motion of Two Objects on Same Trajec- 
tory 

Now, say that there are two objects which satisfy this set of equations. Object, 
1 has the initial conditions which have been discussed, but object 2 trails 
object 1 by a time delay At. Therefore, 

r-i(t) = - At), 

with the subscript referring to the object number. Let the relative motion 
be given by 

Ar = iq - r 2 . 

with the initial state 

r](d)j = r 2 (0 + At) = r 0 . 

First consider the in-plane motion of Equation (11). The in- plane motion 
of objects 1 and 2 is given by 

x t = /1fc.n(Mn)Xu 
x 2 = A k .n(t.-At.J. o)x 0 . 

Therefore, the relative motion is 

Ax = |.'Vu(b /o) — • b. .n(/ — A/, /-o)]x ( ). (1-3) 
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Xow 


Ax 0 = Ax(i 0 ) = x 10 - x 20 

= Af.,n (Uh £<>)xo — Ak,n(Lo — A/.. £o)x () 

= [I ~ AkMt-o ~ At. £ 0 )]x 0 . 

And so, 

xo = [7 — Ak,n(to — A t, £o)] 1 Ax 0 . (14) 

Substituting into Equation (13). 

Ax = to) — Ak,n{t — At, £o)][^ — A k> n(to — At, to)] 1 Axq. 

From the properties of the state transition matrix, 

>U:.n(£ — At, /() ) = Afc,n(£, to)Ak.n(k) — At, to). 

Therefore, 

Ax = [A k ji{t, t 0 ) - A k ,n(t , to)A k: ii(t 0 - At, £ 0 )] 

— • / lfc,n(£o — At, to)] x Axo 

= A k fi(t , to)[/ — Afc.n(tq — At, t 0 )] [/ — Ak.a(to — At, to)] 1 Ax (1 
= A k> n(t, t 0 )Axo- 

The same result may be applied to Equation (12). The resulting relative 
motion equations are then 

Ax = A k 'Si(t, to)Axo 
Az = A w . w (t,t 0 ) Az 0 , 


with the full matrices as given in Equations (9) and (10). 

It may be preferable to present the relative motion in terms of the relative 
position at the time of insertion of the second object — that is, at time to+At. 
Referring to Equation (13), 

Ax(fo + At) = [Afc,ji(fo + At, to) — Ak.ii(t-o, £o)]xo 
= [^r,sj(to + At, to) — /]x o- 
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Following the same logic as above, 

Axo = -A;.: ji(/o + At, to) 1} 1 Ax(to — At). 

and so 

Ax = [A ksl (t : t 0 ) - A kn (t - At, t 0 )]x 0 

= [A k n{t — At, t())A k <i(t 0 + At. to) — A k .n(t — At, to)] 

[^A:.n(/ : o + At, to) — J\ 1 Ax(to + At). 

= A k si(t — At, to)Ax(t 0 + At). 

Once again, the same result applies to the out-of-plane motion, giving 

Az = A^(f - At, t 0 )Az(t u + At). 

Note that these results could also be written as 

Ax = A k <i(t, t 0 + At)Ax(t 0 + At) 

Az = A^ w (t, t 0 + A/;)Az(t 0 + At). 

7.2 Relative Motion With Insertion Position or Time 
Error 

Next, consider the possibility of an injection error for the second object, such 
that it does not in fact have the same initial state as the first. Let Xot refer 
to the in-plane state error at injection, such that. 

x 2 (t 0 + A t) = x n + x 06 . 

The resulting relative state is then, simply, 

Ax = Ah.<i{t, to + Af)(Ax(fo + A t.) a — xo (,), 

where Ax(to + At) a is the nominal relative state at the second object’s in- 
jection time, as discussed above. Again, the same result applies to the out- 
of-plane motion, giving 

Az = to — At)(Az(t 0 • •• A t) a + z 0b ). 
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It may also be useful to examine an error in the injection time of the 
second object. Let the time error be A such that 

A t = A t n + Af/„ 

where A t a refers to the nominal delay between the insertion of the first and 
second objects. In this case, the in-plane state of the second object is 

x 2 = A kM {t - A t a - A t b ). 


The relative state is then 


Ax = [A k si{t, t 0 ) - A kn {t - A t a - A t b , h))l x o- ( 15 ) 

At the nominal injection time, to + A t. a , the relative state is 

Ax(to + A t a ) = [Ak.nito + A t a , to) — Ak : n(to — A tb, fo)]xo- (16) 


The second state transition matrix may be expanded in a Maclaurin series 
about At./, = 0. A first-order expansion gives 


Ak,u(to — A tb, to) 


0 -fA t b 

kQAtb 0 


= I + B k #(At b ). 


Substituting into Equation (16), 

Ax(h) + A t a ) = [A^a(t{) + A t a) to) — i]xo — B^.^(Atf } )x. q. 

Note that the first term gives the nominal relative position at the nominal 
insertion time of the second object; the second term gives the error. Denote 
this first term as Ax(/ 0 + A t a ) a . Then, 

x 0 = [Ak,u{t 0 + Ata.to) - /]“'Ax(fo + A t a ) a . 

Substituting into Equation (15), and again invoking the properties of the 
state transition matrix, 


Ax — [Ak^t, to) — Aksiit - A t a — A tb- fo)] 

[- / U,nCu + A t a ,to) — f] _1 Ax(fo + A t. a ) a 


19 


-- [Ar,u(Mo) — shiSiil- “ A/„, /.p)/!^ u(l() — A/./,, /{))] 

+ A l: a . to) — f] 1 Ax(/ f ) — A/ a ) a 
— {A^oU, to) ~ - U.n(/ _ A/ a . /<))[/ -- Bksi{At b )]} 

\A k M{to + A /: (l , /o) ~ /] Ax(/q H- A/ a ) a 
= Ak,n{t - A t rt) /'o){/l/;,u(/o ■+• A/ a , t 0 ) — / + fifc^A/)]} 

|^A:,S 2 (/-a + A/ f/ , d)) — /| 1 Ax(/() -f A l a ) a 
= A/,: s n(/ — At„, £o){J ~ 5 /^q(A/)[/ 1/, n(/ {) -f A t a , to) — /] 1 } 

Ax(/ 0 + At a ) a 

— At a} <o)Ax(to + At n ) 0 

~ -4/r,u(^ ™ At a , to)Bk^i{Ai) Ak{i(to + A t a . to) — /] 1 

Ax(t 0 + At ri )„. 

Here, note that the first term is again the nominal relative position, while 
the second term is the resulting error. 

8 Some Examples 

Currently, there are four different MATLAB scripts that present numerical 
solutions to these models: 

• full equations of motion about L? 

• linear equations of motion about Lo 

• quadratic equations of motion about 

• relative motion of two objects on the same trajectory. 

We first duplicated Hamilton's results shown in his Section 8.1.1 to check 
the linear code. 

For subsequent analysis, we used a consistent set physical quantities taken 
from Dunham and Muhonen [G]. 

Let us examine the same initial conditions applied to these four numerical 
solutions. We begin near the L-> equilibrium point. Consider displacing 
the spacecraft 10,000 km away from L>z along only the x direction, with 
y(0) = —1. 127.8 km/day as defined from Equation (8) of Section 5. All 
of 7(0), 2 , and 7(0) are zero. As noted earlier, motions in the Z-planc 
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Figure 4: Solutions to Models Show Contrasting Fidelity. 

are uncoupled from the X-Y plane: therefore, for these initial conditions, all 
motions will be in the X-Y plane. 

Figure 4 contrasts the differences among the three models of the equations 
of motion as shown for 170 days (nearly a full period of 178 days). The full 
model does not even close on itself. For these intit.ia.1 conditions the full 
model diverges and the quadratic model shows very good agreement with 
the full model. The linear model does close on itself, but only represents the 
actual motion for about one-quarter of an orbit. Shown on the figure is the 
Li point, which is the origin, and the starting location. MATLAB’s "ode45” 
solver was used, which is based on an explicit Runge-Kutta (4,5) formula. 
The print interval was one day. 

Continuing with the same data produced for Figure 4, we turn to Figure 5. 
Here we take the component positions of the linear and quadratic models and 
subtract the corresponding position of the full model. Although plotted for 
22 days beyond one period, we can see how rapidly the divergence grew as 
the second orbit began. 
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Figure 6: Relative Position of Two Spacecraft. 


We move on to the analytic solution of the relative motion between two 
spacecraft. Here we used the linear solution coded from Equation (13) of 
Section 7.1. See Figure 6, which shows the motion of spacecraft #2 with 
respect to spacecraft #1. The calculations based on the initial conditions 
discussed above were made for each spacecraft. They have the same initial 
position and velocity conditions yet spacecraft #2 began one day after space- 
craft #1. Each diamond on the figure shows the position of spacecraft #2 
in steps of two days. You might consider yourself to be at the origin riding 
on spacecraft #1. With a one-day delay, the position of spacecraft #2 is 
plotted every two days (to offer clarity between plot symbols). At the orbit’s 
beginning and halfway points, the relative change of position of spacecraft 
#2 is less than other orbit points. 
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Appendices 

A Derivation of the Full NonLinear Equa- 
tions of Motion for the Circular Restricted 
Three-Body Problem (relative to barycen- 
ter) 

The development of the equations of motion will follow that of Hamilton. In 
our reading of his thesis to understand his work, we rederived the equations 
of motion. More details are given in these respective appendixes to maintain 
a record of our labors. 

Returning to Figure 1, which is a rotating coordinate system, we find the 
need to add an inertial reference frame to define the rotation vector. 

An inertial coordinate system is defined with origin at the system barycen- 
tcr. The orthogonal basis vectors h\, n 2: and h$ are fixed. Vectors fq and 
a 2 lie in the same plane as the rotating basis vectors d\ and d 2 . The vector 
h : > is always aligned with a 2) . Once every revolution, h { is aligned with a A 
at the same time that n 2 is aligned with a 2 . With respect to the inertial 
coordinate frame, the two large bodies are rotating about their barveenter 
with a constant angular velocity 


0 
0 
U) 

— uJ 71 3 


where 

lG{m i + m 2) 

"W D* S 

i 

and G is the gravitational constant. 

With rotating reference axes, the calculation of the spacecraft's velocity 
with respect to an inertial frame is 


dr 

dt n 


dr 
dt a 


-}“ U) X T 
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where the position vector of the spacecraft in the rotating frame is 

r — A hi ~h Yci '2 Z&s • 

Substitute into the calculation and collect terms: 

a i a 2 a 3 
0 0 cc 

a y z 


t~ — X ci\ -\- Y hi + Za.3 + 

(I t n 


T — (A — 1 i “l - (Y Aa')fl2 “h Z(l%. 


With rotating reference axes, the calculation of the acceleration of vector 
f with respect to an inertial frame is 


cl 2 r 

dt ' 2 n 


- - - dr ct 2 r 

— U) X T a + UJ x (cc X T a ) -K 2 oj x — — H — — — 

ut a (It a 


The first term is zero because the system rotates at constant angular 
speed. 

Substitution into the equation yields 


dt 2 n 


iO x ( — u)\ ci\ y u 2 ATh 2 ) “h 2 (jj x (X"hi -t- V 0.2 4“ Z d-d) 


+ (A hi -r Y(i 2 + Zh:\) 


rtl 

<> 2 

(h 


hi 

«2 

dz 


0 

0 

UJ 

+ 

0 

0 

2 u) 

+ (Aai + Yci‘2 + Za 3 ) 

-ufY 

wX. 

0 


X 

Y 

z 



r = (X - 2 ufY - u! 2 X)d\ + (Y + 2u>X - 1 JY)a, 2 + Za 3 

We associate the kinematics of F with the force of gravity acting on the 
system. Gravity is the only force acting on the system. 

From Newton’s second law, which states that the time rate of change of 
linear momentum is proportional to the force applied. For our fixed-mass 
system of interest 

=, - Gmiin-i _ Gm^m-i _ 

F = m 6 r = : — p - — r 1 : — —r 2 . 

ri 1 r 2 F 
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As a reminder. in \ is the mass of the larger largo body, m 2 is tlic mass of 
the smaller large body, and m : > is the mass of the spacecraft. Now substitute 
into this expression the equations for r, vq and r-> and separate the vector 
equation into the throe components. 


( A — 2 m V — a ; 2 A J o ] f V A 2 mA c m V ) o ) ~\~ Zci% — 

— 7 — -y[(A + Di)ii\ + V a ■> + Z d;j] — " ' . ^[(A — D) )(t\ + 1 a 2 + Za^] 
r, J n J 


A - %jY - m 2 A = - 


Pi (A + D,) p,(A-D,) 


V -l- 2mA - u/V 


2 V' 


/W /nr 


z = - 


|nl :i Ini 3 

PiZ /nA 

Ini 3 In I 3 


From Wie [4] we note these equations of motion can also be expressed in 
terms of a pseudopotential If = U{X, Y, Z) as follows: 


X - 2m Y 


*Lku x 

dX x 


Y + 2mA = % = U Y 
or 

Z = 

()Z 

where f.lio pseudopotontial U, which is the centrifugal plus gravitational force 
potential, is defined as 

U = im 2 (A' 2 + V' 2 ) + — + — . 

2 /'i r 2 

Also 


n = v (A' + D[) 2 + Y- + Z 2 


r-2 = \/(X - D^ + Y'^ + ZY 
Again, the initial conditions to be specified are 

A(0) A(0) 

v( 0 ) y-( 0 ) 

Z(0) 2(0). 
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B Derivation of the Full Nonlinear Equations 
of Motion for the Circular Restricted Three- 
Body Problem Near L 2 

Begin with the results of the derivation of the full nonlinear equations of 
motion relative to barycenter. Equations (1)— (3): 


X - 2 uY - u> 2 X 
Y + 2u>X - lj 2 Y 


Hi{X + D i) {^{X — £> 2 ) 


Z = - 


lor 

Mi Y M2 y 

ln | 3 |n | 3 

Vi z {.12 z 


in 1 


n 


in | 


Substitute using (4)-(6): 


X — Xq + x 
Y = y 0 + y 

z = Z n + z, 


to obtain 


(A'o + x) — 2u(Y 0 + y) — u> 2 (Xq + x) = — 
(>0 + y) — 2 oj(Xq + x) — u 2 (Yq + y) = 


ju i(Aq + x + D\ ) { 12 {Xq + x — D 2 ) 


In I 


Ini 3 


(Zo + i) =. - 


a<i(*o + y) _ /*2(Vo + y) 
In I 3 |n | 3 

{1,\{Zq + z) {I.2(Zq + z) 


n 


n 


By definition the equilibrium point is stationary; therefore, many terms 
are zero. 

X 0 = 0 x = 0 

Vo = 0 y = 0 

Zn = 0 z = 0. 

The equations reduce to the following: 

;; o. ,2 (Y , X _ yi(X 0 + X + Di) {P>(Xq + X - D 2 ) 
x — Zajy — u j + x) — 5 
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y -f 2'jj.v 


— ^ ,2 (Vu + y) 


/ /■ i ( Vo + y ) M2 (Vq + y) 


liii.Za + z) H2(Z 0 + z) 
z ~ 3 r 3 

'l >2 

where ?‘i and r% are now defined below as 

/'] = . \J (A'o + x + D\ ) 2 + (V 0 + y ) 2 — (Zo + *) 
r 2 = y (A"o + — A >) 2 + (Vo + ,(/) 2 + (Zo + ~) 
The initial conditions to be specified are 


.r( 0 ) x( 0 ) 

/;( o) m 



C Derivation of the Linear Equations of Mo- 
tion for the Circular Restricted Three-Body 
Problem Near L2 

We develop linear equations to have models applicable for linear analysis 
methods. 

We document, details not shown in Hamilton’s thesis. 

We repeat the nonlinear equations for convenience. 


where 


Also 


X - 2 loY = 
Y + %)X = 
Z = 


dU 

dX 

du 

dY 

du 

dz 


U = X 2 (X 2 + Y 2 ) + — + —. 

2 r 1 r 2 


r, = ^X + Dtf + Yt + Z* 
r-2 = \f(X - D 2 )> + Y* + ZA 


We substitute the individual coordinates 


X — Aq + 
Y = Y 0 + y 
Z = Z() + z 


to obtain 


(A'o + x) — 2uj{Y 0 + y) 
(Vo + y) “ 2lu(Xq + x) 


(Z 0 + z) 


°± = v x 

oX 

■ou _ 

dY - Uy 

oz z 
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Xow, perform a Taylor series expansion on these parti als of U, about the 
libration, equilibrium, point. The first two terms of this series for a function 
of one variable are 

/{■*)& f(X [ >)+.vf'(X t} ). 

For U x , \et f(x.y.z)=^ 

Note that because the function is already a first derivative, we stop the 
Taylor series at its first derivative so we can stop the series at the quadratic 
term. We have a function of three variables so the Taylor series is modified 
as follows (with L 2 substituting for subscript 0 as the chosen expansion and 
evaluation point) 

f(x,y. z) = U x | 

For U y , let /(; r. y, ~) = ^4, then 
f( X ,y :Z )=U Y \ Li + X^Uy 

For U z , let /(. r, y, c) = then 
f(x.y,z)=U z \ L2 + x~U z 

Bv definition at the equilibrium point, there is no motion. Many terms 
are zero. 


l 2 


d 

+ VgyUx 


+ z—Ux- 
l> 9Z X L2 


L-2 


+ Vgy' UV 


Li 


+ z m Uy 


0 () 
+ !J W bz + 

L 2 L-2 




Y i. 2 = .i'L 3 = U x \ l-2 = 0 
U, = ()l 2 - Uy\l 2 = 0 
Zl> = z l 2 = Uz | l 2 = 0 

Additional ly. terms of U\ and ! V involving t.ho partial of Z are zero be- 
cause motions in the Z-direetion are uncoupled from the X-Y plane motions. 
Therefore 


x - 2u>y = x(J x .y|l 2 + yU\x\L 2 
y - ‘hoi: = x0.xy\l 2 + vUyy\l 2 
•' = z Uzz\l 2 
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The above equations are the first-order, or linear, coupled, equations of 
motion for the circular restricted three-body problem. These describe the 
motion of the spacecraft relative to an equilibrium point (Xo,Y( h Zq). The 
equilibrium point rotates with a constant angular velocity. 



For the restricted three-body problem at £2, the Yq and Zq terms are zero 
and the equations reduce to 

z-2ujy = xU X x\l 2 
y + 2u;i = ijUyy\l 2 
z = zUzz\l 2 
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Return to the general form of U as we first carry out the derivatives of 
U, then substitute for the libration point. 


U = ~(X 2 + Y 2 ) + 


{. x + A) 2 + y 2 + z 2 j | J(x - Do y + y 2 + zA 


For X: 


8U BU IV 2 J 

•Tr 7 = ^ ) 


dX OX 2 


X + A) 2 + Y 2 + Z 2 


0X I \/(V - D 2 ) 2 + Y 2 + Z\ 


+ / 1| ||[((A- + D l ) 2 4-r 2 -Z 2 )^] 

+ B2^[((X - D,) 2 + Y 2 + Z 2 )-^} 

Employ the differential d(u") = uu n ~ l du to calculate the last, two terms 
of the equation : 


( W + Ci) 2 + V 2 + Z ? ) 2 


/'i( — o)(( y ^ + Ci) 2 + V ’ + Z 2 ) 2 — ^r((X + A)" — c 2 _ Z") 


= , =- (2X + 2A) 

2 S J(X + A) 2 + V' 2 + Z 2 p 

/'•iA + A) 

^(x + A) 2 + V' 2 + x 2 ) 3 

Of course this calculation is similar for the term beginning with ju*. Here 


assign 


X + A ) 2 + V' 2 + Z' 2 


X - Do) 2 f V 2 -f Z 2 
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Finally 


du 2v m {X + Di) MX-d 2 ) 

— uj A — — 


flA 


rf 


ro 3 


Now take the second derivative using cl 


udv + vdu: 


+ ((X + D,) 2 + r 2 + Z 2 )-i)] 


d 


<9 


— ^[(X - D 2 ) ((A" - D 2 ) 2 + Y 2 + A 2 ) *)]. 


^ = u; 2 -fi 1 i [ (X~D 1 )\ i -~((X + D l f + Y 2 + Z 2 )-i(2X + 2D l ) 


+ ((A + D,) 2 + r 2 + Z 2 )-5(l) 


1^— ^((A' — D 2 ) 2 + V 2 + A") - (2 A' — 2 Do) 


- & | (A - Do) 

+ ((A - /; 2 ) 2 + f 2 + z 2 )-5(i)J. 
Now reduce the terms and substitute for rj and r 2 


0 

'1 3(A + D,) 2 


1 3(A - D 2 ) 2 ' 

LlT - fll 

r 3 r 5 

L ; l 'l J 

- l '-2 

r\ r 2 


For Y: 


BU ., v . //iF //. 2 F 

a? = “'' “ 7f 


and 


■> 

' 1 3F 2 


' 1 3V' 2 

a;" - /il 

r 3 r 5 

L r i r i J 

CM 

-< 

! 

r '2 r\ 


For A: 


9D 

iJA 


fiZ fi 2 Z 


'i 


r I 


and 


Uzz = —l>i 


' 1 3 A 2 


' 1 3 A 2 

3 ...5 

- /<2 

„3 ,,5 

l r i J 


L r 2 } 2 J 
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D Derivation of Quadratic Differential Equa- 
tions for Circular Restricted Motion Near 
L 2 


We repeat the nonlinear equations for convenience: 


where 


and 


.V - 2 uY = 
V + ‘2uX — 
Z = 


du 

dx 

du 

dY 

du 

dz 



(17) 

(18) 

(19) 

( 20 ) 


n 2 = (X + D x ) 2 + Y 2 + Z 2 (21) 

r 2 2 = (X - D 2 ) 2 + Y 2 + Z 2 (22) 

The right sides of the equations of motion can be formed from differenti- 
ation of Equation (20): 


9U_ _ 2 Hi dr i /tg dr 2 

dX Ty 2 dX r 2 2 dX 

dU 2v . A*i dry A 4 2 dr 2 

dY = u ~ 'rfdY 

8U Hi dr i // ; Or 2 

dz = ~V{-~dz ~ ^-~dz' 


(23) 

(24) 

(25) 


The required derivatives of rj and r 2 are formed by differentiating Equa- 
tions (21) and (22): 


dr x _ X + D i 

dX ~ Vy 

On _ y 

ay “ 7y 

dry _ z 

JZ ~ Vy 
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Equations (23) - (25) are expanded in Taylor series expansions about the 
L 2 libration point, at coordinates (X ,Y,Z) = (A’o, 0.0). Let 


V 


X - X 0 
Y 

Z 


(26) 


give the displacement from £ 2 - Denote DU /OX by (J.x ■ and similarly for Uy 
and Vz- Then, the Taylor series through second-order in x.y,z are given by 


Uy = U 


.Vi L, 


du 


X 


dx 


X + 


OU x 


c)Y 


y + 


DU \ 


DZ 


1 ( D'HIy 


d 2 U x 


D 2 U X 



A vertical bar with subscript. L 2 indicates evaluation at the coordinates of 
the libra, tion point. 

The derivatives of U are constructed by differentiating Equations (23) - 
(25). Note that higher-older derivatives of rq and r 2 arc required. These are 
given by 

d 2 n 1 (X + DO 2 ' 

dx 2 ~ n n 3 

dY 2 rq rq 3 

d 2 rq _ 1 Z 2 

Jz 2 ~ 

d 2 n {X + Di)Y 


dXdY 


ri 3 

d 2 n 

(X + Di)Z 

dXdZ 


n 3 

d 2 r 1 


dzoz 

7*1 3 


d 2 r-. 

1 

(X - £>2 

dx 2 

Vo 

r 2 3 

d 2 r 2 

1 

F 2 

dY 2 

r 2 

T2 3 

d 2 r 2 

1 

z 2 

dz 2 

? 2 

r 2 3 

d 2 r-2 


- D 2 )F 

dX dY 

= 

?’ 2 3 

d 2 r 2 

_(* 

-£ 2 )Z 

dXdZ 


r 2 3 


8 2 r 2 _ YZ 

dZdZ ~ r 2 3 

The evaluated partial derviatives are then: 


dU x 

2 1 

2/X2 

dX L 

W (yY 0 + £i) 3 

(* 0 - D 2 ) 3 

dU Y 

. .2 /h 

/<2 

dY l 2 

(Xq + £| ) 3 

(Z„ - E> 2 ) 3 
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and 


dU x 

Hy 


du z 



/'i 

l lt 2 

dz 

L, 

(X„ + A) 3 

(A 0 — D?) 

d 2 u x 



6/4 

6/4 

ax 1 

1* 2 

{x 0 + D 1 r 

(X 0 - D 2 ) 

a 2 u Y 



3/4 

3/4 

i)Y 2 

L> 

(X 0 +/M 4 

(X 0 - lhY 

d 2 U x 



3/4 

3/4 

az 2 

L>2 

(Xq + d,y 

(X 0 - D 2 ) 4 

dUx 


dU Y 

_ iPUx 

a 2 f/.v 

DZ 

Lj 

DZ 

, ~ OXY 

L o 

“• axz 


a 2 u x 

d 2 u y 

(f-Uy 

0 2 U y 

d 2 U z 

avz 

. ~ av 2 

L-2 

, ~ dYZ 

L-j 

X 822 

l 2 ~ 


= 0. 


Note that, because of the continuity of the derivatives, any other derivatives 
are, in fact repeats of those listed. For example, 


d 2 Uy 

a 2 ( 

r au\ 

ax 2 

ax 2 \ 

K dY) 


aHi 



ax 2 Y 

a 2 

( au ' 


axv 

[ax 


a 2 tix 



QXY 


By definition, the L -2 coordinates satisfy the equilibrium condition of 
Equations (17) - (19). Therefore, 

tfvL = u y \ Li = u x \ b2 = 0 . 

Substituting the derivatives into Equations (27) (29), 


Ux = 


+ 




+ 


t/12 


(A 0 + D\Y (X Q - D 2 f 
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U Y = 


U z = 


Define 


3 

Mi 

M 2 

J 

~~ 2 

[(Xo + D^ 

(X 0 -D 2 ) 3 J 

, ,2 

Mi 

/*2 


Lis ~ 

(X 0 + !h ) 

(X 0 - D 2 f J 

+ 3 

Mi 

, / 7 '2 

J 

_(x 0 + D 1 Y 

(X 0 - d 2 ) 3 J 

\ 

Mi 

M2 

2 

(* 0 + A ) 3 

(A'o - Do) 3 

+ 3 

Mi 

/^2 

J 

_(X„ + A) 4 

(X 0 - D 2 ) 3 J 


(2.r 2 - y 2 - ^ 2 ) 


y 

xy 


xz. 


A k 


B = 


Mi 


(X 0 + A ) 3 + {X Q - Do ) 3 

Mi M2 


M2 


(Xo + DiY (X 0 — D 2 ) 1 
Then, substituting in Equations (30) - (32), 


Dx = (^ 2 + 2,4).r - ±B(2x 2 - y 2 - z 2 ) 

Uy = (uj 2 - A)y + 3Bxy 
Uz = -Az+Wxz. 


(30) 

(31) 


(32) 


(33) 

(34) 


Substituting into Equations (17) - (18), and using Equation (26) on the left 
side, gives 


x — 2uiij - or x — 2 Ax — ^B(2x 2 — y 2 — z 2 )- 

y + 2ux — 'J 2 y — —Ay + 3 Bxy 
z = —Az 4- 3 Bxz. 
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Part 2: Linear and Quadratic 
Modelling and Solution of the 
Relative Motion 
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1 Introduction 


Some continuing long-term goals of our sponsor are to develop high-fidelity 
equations of motion representing the formation flying of spacecraft near the 
Sun-Earth L 2 point and equations describing the relative motion between 
these constellation members. The equations are to be used to develop orbit 
control schemes. 

This is our second report, which continues from Part 1 [1], with further 
derivations of analytical expressions for the relative motion between a for- 
mation of spacecraft orbiting near the Sun- Earth L 2 point. This research is 
loosely motivated by formation flying concepts for the MicroArcsecond X-ray 
Imaging Mission (MAXIM) Pathfinder. 

To begin the understanding of the behavior of objects in formation near 
L 2) we begin with the assumptions of a circular restricted three-body prob- 
lem. In this report, we are not modeling, the true eccentric orbit of the Earth 
and Moon about their respective primaries and other planetary and solar 
gravitational orbit perturbations. 

In this report, we present a preliminary understanding of the relative 
motion between two spacecraft orbiting the chosen L 2 point. While com- 
ponents are further identified below, we describe the motion of a typical 
telescope spacecraft with respect to the central hub spacecraft of the constel- 
lation. The hub is treated as the constellation’s reference point; its orbital 
path about L 2 was examined in Part 1. The description of motion of one 
telescope spacecraft with respect to the hub spacecraft can be applied to as 
many telescope spacecraft as needed. 

After a significant literature search, we state that what we are uniquely 
contributing is a description of spacecraft in formation flight about the Sun- 
Earth L 2 point. 

First, we list the components of the MAXIM mission. Second, we develop 
the differential equations of motion for the telescope spacecraft with respect 
to the hub. The resulting solution is developed using the Lindstedt-Poincare 
perturbation method to ensure periodic motion. Periodic motion of the in- 
plane and out-of-plane orbit is desirable to maintain the formation. Finally, 
we provide analytic solutions of relative motion. This report thus contains 
the following: 

• development of differential equations of motion for the telescope with 
respect to the hub 
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• description of the linear telescope motion about the hub 

• discussion of the preference of a halo orbit (rather than a Lissajous 
orbit) of the hub about the L 2 point, 

• inclusion of linear hub motion effects upon the motion of the telescope 
relative to the hub 

• inclusion of quadratic hub motion effects upon the relationships be- 
tween the telescope frequency and amplitude 


2 Problem Definition 

The current configuration of the MAXIM Pathfinder mission is depicted and 
described below. 



r 2 


\ 


r 3 


Figure 1: Aperture Formation of MAXIM Pathfinder, Concept of June 2002. 




Seven spacecraft shown in Figure 1 are in the same plane forming a flat 
aperture. The optical hub is in the “center” of six free-flying spacecraft 
(telescopes). Dashed lines do not indicate a physical connection, but rather 
indicate random radial directions listed rq through r 0 . The scalar length of 
r ranges from 100 to 500 meters. Additionally, the six free-flying spacecraft 
are loosely spaced 60° from one another. 

There is also a detector spacecraft located approximately 20,000 km and 
90° out of the plane formed by the flat aperture. The entire system is in 
a Sun-Earth L 2 orbit, which has yet to be designed. The analysis of the 
detector’s orbit relative to the aperture plane was beyond the scope of this 
work. 

At the beginning of this task, we were directed by the sponsor to inves- 
tigate the following: 

1. types and fidelity of models used to describe the relative orbits of the 
seven spacecraft forming the aperture 

2. relative orbits of the seven spacecraft forming the flat aperture 

3. slewing motions of the flat aperture 

The optics hub slews only in attitude. The entire system points all over 
the sky — there is no nominal inertial orientation. Also, we may select any 
distance between L 2 and the optics hub. Additionally, the question was raised 
as to whether the motions of the telescope relative to the hub be described 
in cylindrical coordinates. This may be easier to describe when working with 
the project astronomers. 


We were unable to do all of this with the funding available for this phase 
of the project. The work presented here covers the depth needed for items 1 
and 2. We would need to continue to extend the development for one tele- 
scope into that for six telescopes. A computer-generated visualization of the 
aperture plane’s motion based on our equations of motion would be highly 
useful. 

We have concerns about being able to arbitrarily select the distance' be- 
tween L >2 and the optics hub. We suggest that the MAXIM Pathfinder mis- 
sion may best work in a halo-type orbit. A halo orbit is defined when the 
frequency ratio of out-of-plane to in-plane periods is a rational number. The 
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three-dimensional halo orbit doses on itself and is periodic, in contrast to 
a, Lissajous orbit in which the trajectory does not close. In particular, we 
recommend a halo orbit in which the in-plane and out-of-plane periods are 
equal. 


In Section o we form two sets of differential equations of relative motion: 
full nonlinear and truncated nonlinear equations. 

In Section 4 we present an analytical solution to these equations. 

In Section 5 we show simulations. 

In Section 6 we summarize this report. 


3 Relative Motion Differential Equations 


In Part 1 of this research [1], the general second order differential equations ol 
motion were constructed for an object near the Sun- Earth L-> libration point, 
using the force model of the classical circular restricted three body problem. 
In this model, the Earth is treated as being in a circular orbit about the sun, 
the spacecraft mass is considered to be negligible as compared to the two 
primaries, and only point-mass gravitational forces are considered. 

For this system, depicted in Figure 2. the differential equations of motion 
for an object (object /) near the Sun-Earth Lo are given by 


r i = 


('ilL + Jllj r , _ ( P'l + ^ - M ) x + 

V/'H 3 l>2 i A J V l> P2, ) 


where 


r ( = vector from L 2 to object i 
//■! = solar Keplerian constant 

/(, = terrestrial Keplerian constant (Earth + Moon) 
p u = distance from Sun to object i 
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Figure 2 : Coordinate Axis Definition. 


p2i = distance from Earth-Moon barycenter to object i 
= distance from system barycenter to L 2 
D x — distance from system barycenter to sun 
D2 = distance from system barycenter to earth-moon barycenter 
x = unit vector parallel to sun-earth line of syzvgy, 
pointing in sun-to-earth direction 
n = terrestrial mean motion about sun (assumed constant). 

Let r h and r t denote, respectively, the vector from L 2 to the hub and to 
a telescope. Therefore, if r is the vector from the hub to the telescope, the 
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differential equation of motion for the teleseope relative to the hub is 


r = r t - r,, 
AM 


Pit 3 Pat 3 
Pi | /'2 


rt _ | Pl( ;f: e + ^l) + l>2(,Ve - Ai) f x 


Pit' 


Pat " 1 


3 


= “Pi 


Pi/, P 2/1 

rt r h 


r/> + 


/ti(;r e + Di) P 2 (x e 


D, 


Pit 


3 


pi/P 


/'■a 


Pi li. 

Yt 

pat 3 


3 


P 2/1 


T/t 

P2/1 


Pl('£'e + -Dl) 


1 


i 


Pi / 3 Pi /, 3 


x - fi 2 (x e ~ D 2 ) { 


1 


X. 


( 1 ) 


v AMb Plh 

where now, in general, the subscripts h and t refer to the hub and telescope. 


3.1 Series Expansion of Differential Equations 

The differential equations of motion are now expanded in powers of the dis- 
tances between L 2 and the hub, as well as between the hub and the telescope. 
This is done with the intention of developing an analytical solution in terms 
of these quantities, useful for performing a control system analysis. 

First, the inverse powers of the p magnitudes are written in terms of these 
distances. If p lh is the vector from the sun to the hub, then 

P\h = ( ' : e ~ D i)x - r h . 


The square of its magnitude is then given by 


P\h 


2 


Plh ’ Plh 

(:T e + D\Y -b I'k 2 — 2(x e + D\)( r/ 7 • x) 
(;r e -f D\Y + r/P Hb 2(.r c -b D[)xk 


(x r + /)] ) 2 



2 :r, 

x e + D\ 


where X;, is the x-component, of r /,. Them 


1 

PvY 


(x e + D \ 


1 + 


'i'h 


X e -p D\ 


+ 


2 xk 
x e + D\ 


i - 3/2 


— Cri T & 1 ) ’*(1 ‘I c \h) 
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whore 


*1/1 


a / r h 


+ 


2 x h 


x e -f D\ J x e + D\ 
is assumed to be less than unity. Using a. binomial expansion, 

1 1 * 

Pih 3 ( x c + D t ) 3 X o V k J lh 


1 


(x e + D x y 
In the same fashion, for the telescope, 


-t(T) 

k - 1 v 7 


(Mi 


1 


1 


Pit 3 i x e + -Dl ) 3 


i+Elfk 1 


where 


A 

— 




A:— 1 


+ 


2x t 


x e + Z)i / x e + Z?i 

again assumed to be less than unity. 

Xow, r t may be written in terms of ?> t . Using 

r t = r h + r, 


( 2 ) 


Also, 

Therefore, 


n 2 = (^h + r) ■ (r/, + r) 
= r h 2 + r 2 + 2r h • r. 


Xt = X h + X. 


eu 


n 


+ 


2.r t 


Xe + D\ J X e + D\ 

'i'h. 2 + t 2 + 2r/, • r 2(.c/, + x) 


+ 


(x e + D i)’ x e + D\ 


i\ 


+ 


2x h 


(Xe + Di) 2 X e + D\ 

= cj h. — Si, 


+ 


r 2 + 2r/, • r 2x 
_ (x e + D\) 2 x e + Di 
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where 


^ a f + 2r /t • r 2x 
(;r f . -f D]) 2 x r + 0\ 


Again using a binomial expansion, the powers of e\ t , as required in Equa- 
tion (2), are given by 

tu — Ui/i + &i) k 

k 


E 

f :- 0 


a A:-* 
r j / i 


1/1 


+ E 


\X A A;-* 
01 e l/r 


This expression may now be used in Equation (2). giving 


1 


.1 


Pi/ 3 (x t: + D { )* 
1 


1 + 



A— L 


) 3 

1 


i + E 


A.-1 


-3/2 

A: 

-3/2 

A: 


^l/l 


E 


\ a £ . k-r 

a Ml e l h 


f-'l/i 


' ' l; A-=l x 7 C=1 


pECf)E 


I r /!, k—t 
Ol fl h 


1 


»E 


-3/2 


Pi/! 3 (x t + Di) 3 ^V k 

Therefore, for use in Equation (1). 


i:= 1 


E 


Pit 

Additionally. 


3 


Pi /,/ 3 (Xe + A) :i f^ V A: 


e(T)ec 


l a- C k-t 

Ol c \h 


f=l 
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r/,. r h + r r h 
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Pi/"’ Pi/T 
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Pi/ J 


( 4 ) 




An identical development may be used to form the following relationships 
for distances involving the earth: 


and 


1 

P2t 3 


1 

P'2/i 3 


1 

{x e - D 2 ) 


E 

k = 1 


-3/2 

k 


E 

f=i 



ft 

P2t 3 


P2h 3 


= r h 



1 N r 

P2/1 3 / + P2t 3 ’ 


where 


(5) 


(6) 


A 

t2h ~ 



2 

+ 


2x h 

X c — D'l 


and 


a r 2 + 2r h , • r 2x 

02 = ?TW + 


(x e - D 2 ) 2 x e - D 2 

Substituting Equations (3) (6) into Equation (1), the telescope motion 
relative to the hub is then given by 


r = — //. t 
- M 2 
= “Mi 


(r ,. + («,+ »,)*) 

(r ‘ + (t * - *>*> G? " E?) + K?\ 


r h + (-^e + D\)x /— 3/2 


(:r e + Z?i) ; 


FE 


fc=l 
oo 


+ 


sE 


it 


-3/2 


'k 


E 


r=i 


- M2 


(x e + Di) 3 ^rV k 

r h + (x £ - D 2 )x 

k = 1 


E 






!h? 


E 


-3/2 

A: 


V / Mx E a.—* 
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(x e - D 2 f 


k=0 


ECflE 


fc 


| ^2 ( j2h k ( 


C = 0 
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3.2 Magnitude Ordering 


A magnitude ordering system is now employed in order to truncate Equa- 
tion (7), prior to forming the solution. It is noted that the motion of the 
system takes place on two separate distance scales. There is a large distance 
scale, in which the motion of the hub about L 2 is described, relative to the 
motion of L 2 within the context of the three-body problem. Then, there is a 
small distance scale, in which the motion of the telescope about the hub is 
described, relative to the motion of the hub about L 2 . 


Table 1: Distance Ratios and Basic Accelerations 


r/r h 

8.333 x 10“ 7 

r h/ ( x c + ) 

3.971 x 1U -3 

7',/(.T e - Do) 

3.980 x lO" 1 

r/( x e + Di) 

3.309 x 10-° 

r/ {x e - D-2 ) 

3.316 x 10“ 7 

l l \/{ x e + D\) 2 

5.812 x 10~ 6 km/s 2 

ni/{ x ( - D 2 y- 

1.775 x 10“ 7 km/s 2 


For ordering purposes, consider the hub motion about L 2 to be on the 
order of 600,000 km, and the telescope motion about the hub to be on the 
order of 500 m. Using the constants of Appendix A, the relative distances 
are approximated by the ratios of Table 1. In the differential equations, 
these ratios scale the basic acceleration quantities which also appear in the 
table. Terms involving these basic accelerations scaled by 77 J{x c + D\) and 
r'h/{ x c ~ D 2 ) are now designated as being of order 1 : terms involving the 
basic accelerations scaled by r/(x e + EV) and r/(x e — Do) are designated as 
being of order 3. 

Retaining terms in Equation (7) only through order 3, substantial algebra 
gives the truncated differential equations as 

f = A [— r + 3xx] 

+ B [3;rr>, T 3x h v + (3r /( • r - 15:c.c/,)x] 

+ C [(3r„ • r - 15.X.T, K + 2 ( Vh 2 - 5x, 2 )r (8) 

- ^(2x h xu • r - 7 xx h 2 4- xt, 2 )x] , 
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where 


Ri 

+ 


M2 

0) 

( x e + D\ ) 3 

{X e 

- d 2 ) 2 

Mi 

+ 


M 2 

(10) 

(x e + D \) 4 

(Xe 

- D 2 y 

Mi 

+ 


M2 

(11) 

(;r e + Di) 5 

(X e 

- w 


Note that this truncation includes terms which are linear in the coordinates 
of r and no more than cubic in the coordinates of r Terms involving A . 
R, and C are, respectively, of orders 3, 4, and 5; lower order terms do not 
appear. 

The acceleration vector r may be written relative to a rotating coordinate 
system which rotates at the constant angular rate n about the z-axis normal 
to the ecliptic, and with the x direction as previously defined. This gives 



x — 2 m) — n 2 x 
ij + 2 nx — rry 

z 


where the column vector notation is used to indicate the xyz vector compo- 
nents. 


3.3 Linear Telescope Motion about the Hub 

Consider now the linear motion of the telescope about the hub. From Equa- 
tion (8), the linear differential equations are given by 

r = A [— r + 3rx] , 


or, in component form, 


x — 2 ny - (n 2 + 2A)x 
ij + 2 nx - (n 2 — A)y 
z + Az 


= 0 . 


(12) 


As would be expected, these differential equations take exactly the same 
form as the linearized hub equations discussed in [1], with the same funda- 
mental frequencies; the same solution approach may be followed. Clearly, the 
out-of-plane motion in the z direction is decoupled from the motion in the 
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x\j piano. The z equation describes simple harmonic, motion, with a solution 
that may be written as 

2 = A z si n(ut + ip), 

where v 1 — .4. 

Following the standard linear analysis, the in-plane motion is assumed to 
exhibit a solution of the form 



Accordingly, the fourth-degree characteristic equation for s is 

s 4 - (A - 2 n 2 )s 2 - (u 2 - 2 A)(A - n 2 ) = 0. (13) 

Keeping in mind the relative magnitudes of n (ss 0.0172 racl/day) and A (~ 
0.0012 rad /day) , the solution to this equation has four distinct roots. One 
root is positive real, and corresponds to a divergent mode: 



A second, negative real, root corresponds to a convergent mode: 

Sc 

The remaining two roots are a purely imaginary conjugate pair, correspond- 
ing to oscillatory motion with natural frequency A given by 

A 2 = - 4 + n 2 + y^(4 “ ft 2 ) 2 + (ft 2 T 2 A) (A — n 2 ). 

Each mode shape is described by the eigenvector associated with the corre- 
sponding eigenvalue s. 

As in the analysis of [1], periodic motion occurs if the in-plane initial 
conditions are selected so as to excite only the oscillatory linear modes. From 
the solution of the eigenvalue problem, placing this requirement- upon the 
initial conditions gives 

• i: 0) = pj( °) 

//(O') = —kXx(O). 


G2 



whom 


, A 2 + n 2 + 2A 
k = — • 

Under this relationship among the initial conditions, the complete linear 
solution for the telescope motion about, the hub may now be written in the 
form 


x = —A x cos (A t + ({> ) 
y — kA x sin(A t + <f>) 
z = /l-sinfW, + VO* 


3.4 Halo Telescope Motion 

It is desired that the constellation of telescope spacecraft remain in approx- 
imately the same specified plane over the few days 7 duration required for 
observations. To achieve this orientation, a halo-type orbit of the telescopes 
about the hub is selected. Such an orbit provides periodic motion in the 
aperture plane, with the out-of-ecliptic-plane fundamental frequency equal 
to A, the fundamental frequency of the xy-planar motion. 

In order to do this, the inclusion of higher-order forces is used to ad- 
just the out-of-plane fundamental frequency. Consider the ^-component of 
Equation (8): 

z=:~Az 

-f B(‘Axzk + izx-h) 

+ C(-I2xx h z h + 3yy h z h - 6zx h 1 + \zy h 2 + \zz h 1 ). 

Recall that A = i/ 2 , and define A such that 

A = A 2 - v 2 . 

Then, this differential equation may be written as 
z + A 2 z = A z + B{oxZ), + 3 zxh) 

+ C(—l2xXhZh + 3yyi,Zh — Qzxh" + + § zzt 2 )- 

Hero, the magnitude of A allows the term A z to be treated as a higher 
order term, grouped with the terms containing the coefficient C. Using this 
formulation, the linear solution now becomes 

: = .b sin(A/ + ip), 
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with the same fundamental period as the in-plane morion. 

Together with the x and y components of Equation (8), the system dif- 
ferential equations are 

x — 2 ny — n 2 x = — 2 Ax 

+ f3(—Gx:i-k + 3 ytjh + 3.3 z^) (14a) 

+ <T '(12.x r/,' 2 - 12 yxhVi,. - I2zx h z h - ii.r.y h 2 - 6 xz h 2 ) 
y + 2nx — n 2 y — — Ay 

+ B(3xijk + 3 yxh) ( 141>) 

+ C{-I2xx h y h + hjyh 2 + 3zy h z h ~ 6 yx h 2 + yz h 2 ) 
z + A 2 ^ = Az + B(3xzu + 3zx h ) ^ 

+ C{-\2xx h z h + 3yyi,zi, — 6 zxt, 2 + 5 z Vh + f- 2 // 2 ) • 

4 Lindstedt-Poincare Development 

The analytical solution to the expanded equations of motion of Equations (14) 
is now developed using a modified version of the Lindstedt-Poincare method. 
This method provides for the sequential solution of a system of differential 
equations, ordered by magnitude of the terms, while simultaneously placing 
restrictions on the initial conditions, in order to ensure periodic motion. In 
this development, the equations are expanded through third order in the 
small quantities, which is defined as terms that are at most linear in the 
motion of the telescope relative to the hub. and quadratic in the motion of 
the hub relative to L 2 . 

Introduce the non-linear frequency terms by scries expansion, and change 
the independent variable from t to r, where 

T - Ojl. 

and 

U.- — 1 ~b uJ\ -f- jJo -(-••• , 

assumed to be an asymptotic series, is used to scale the linear frequency. 
Using primes to denote differentiation with respect to r, the left side of 
Equations (14) become 

J 2 x” - 2 tuny' - n 2 x 
ury" -f 2 unx f — n 2 y 
uS z z" + X 2 2 
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Recall that the forcing terms in the differential Equations (8) are ordered 
such that the terms involving A , B , and C are, respectively, of orders 3, 4, 
and 5. Now, for notations! convenience, these terms are reordered as orders 1, 
2, and 3. The solution vector is assumed to take the form of an asymptotic 
series, such that 

X\ + X2 + r*3 + * * • 

IJl + U2 + 2/3 + ' * ' 1 

+ *2 + Z'3 H 

where the ordering by subscript is consistent with the reordering of the terms 
of the differential equations; the order of a given term is specified by the 
subscript. 

This expansion is then substituted into the differential equations of Equa- 
tions (14), and terms collected by order. The resulting first-order equations 
for xq, jq. and z ly are given by 

— 2 Y7,v/ j - (n 2 — 2/1 ):/■] 

y'l - 2 nx\ + (.4 - n 2 )yi = 0. (15) 

z'[ + X 2 Z 2 

Note that these equations are identical in form to the linear equations of 
Equation (12), with A replacing v in the 2 component equation. 

The second order equations contain contributions from the motion of the 
hub about L 2 . As shown by Richardson [2], the hub motion relative to L 2 
may also be expressed in a series form: 

Xh X\h + X'2h 4 - • * * 

1 !h — D\h + V'lh + ■ • ■ 

Z h _ [ Zlk + 2 2 /i + * * ’ 

In the second order equations, only xu u yih, and z\h are included; at third 
order, the 2 h terms contribute as well. 

The second order equations for x 2 , 7 / 2 , and containing terms which arc 
linear in the position of the hub, are then 

A - 2 ny! 2 - (n 2 + 2A):i: 2 
y” + 2nx',, + (.4 - n 2 )y 2 
4 + A 2 z 2 

-2wi(.e" - ny[) - 6Bxix Vl + 3B{y l y xh h + z x z ih ) 

—2uj\(y” + nx\) + 3B(xiyih + y\Xih) ■ (16) 

— 2 U\z’{ 4 - 3B(x\Z\h + z\Xih) 
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Horn, the telescope motion terms on the right, side are now assumed to he 
known from the solution to the linear equations. 

Next. the third order equations for x->, 1/3, and z%. containing teims lineal 
in the 2h hub position terms and quadratic in the 1//. terms, are 

- 2 ny' s - (n 2 + 2/l).r 3 = 

- 2tx l (x' 2 - m/ 2 ) - 2u^(x'( - ny \ ) - u}\x'[ 

+ 3B(y 2 yih + z 2 z ih ~ 2.T2*ih - 2.ri x 2 h + yym + = 1 - 2 h) 

+ €>C[x i(2x 2 h - y] h - z\ h ) - 2yiXikVih ~ 2z x x\hZ\h] 

2 nx' :i + (A - n 2 )y :i = 

- 2ui(H2 + nx 2) 

+ 'iBlXolJlh + ]j2 x lh + x l lj2li + y\ x 2h ) 

+ 3 C [— 4;ci'C'i/ t ,y]/i + y\ (~2xf h + *y xh + k z (h) + z iVth z ih] 


(17a) 


// 

>h 


2 ^ 2 ( 1 /] + nx \ ) — ^T//i 


(17b) 


J> , \’2 „ _ 

<-3 + A ^3 — 


(17c) 


— 2 uj\ z<? — (2ll?2 + w\)z{ + A27 
-+■ 3 B(X‘)Z\h "b 2*2 X[h + + Z\X‘2h) 

+ 36’ [-4.C|X|fcSi/ t + V\V\ hZ\h + (-2*4 + y\ h + ■ 

Again, the right side terms are assumed to be known from the lower older 
solutions. 

Note that the left side terms at each order take an identical form. There- 
fore, a particular solution may formed for a general periodic forcing function. 
That solution may then be applied to each of the actual forcing terms, and 
the complete solution given by superposition. Higher order homogeneous 
solutions are not required, as they are identical in form to the terms of the 
linear solution, and thus arc already present in the complete solution. 
Accordingly, consider the general case of the forced system, given by 

(18) 

This vector is representative of all forcing terms which occur, where q and 
p arc integers. Next, following the method of undetermined coefficients, the 
particular solution is assumed t.o take the fonn 


' x" - 2 ny' - (n 2 + 2 A)x 1 


Acosiqr + 6) 

y" — 2 nx' - (A - n 2 )y 

= 

B: sinfqr + $) 

z" + X 2 -. \ 


T sin(pr + 7 ) 


X 


R x cos ((/7 + 9) 

y 

= 

R y si n(qr + 9) 

z 


R z sinfjrr + 7 ) 
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Substituting into tho general system of Equation (18), and solving for the 
solution amplitudes, 

2 nqB — A(q 2 + n 2 — A) 
x q A + (.4 — 2 n 2 )q 2 — (A — n 2 )(n 2 + 2A) 

= 2nqA - g(g 2 + n 2 + 2/1) 
y <y 4 + {A - '2n 2 )q 2 -{A- n 2 )(n 2 + 2/1) 

Note that the denominator of R x and R y is, in fact, the characteristic equa- 
tion of the linear system, as seen in Equation (13), where s 2 has been replaced 
by —q 1 . Therefore, upon examination of this particular solution, either q or p 
equaling A corresponds to the resonance condition. In that case, the forcing 
function has the same frequency as the homogeneous solution, and so the 
particular solution must instead have an amplitude which is secular in r. 
However, in order for the solution to be bounded for all r, such secular terms 
may not appear. Therefore, restrictions are placed on the solution, such that 
secular terms do not arise. 

For the case that p = A, it is required that T — 0. For q = A, the 
general x and y-component differential equations are combined into a single 
fourth-order differential equation: 

x"" - (A - 2n 2 ).x" - (A - » 2 )(n 2 + 2.4 ).r = 

[2 nBq + (.4 - n 2 — q 2 )A] cos (qr + 0). (20) 

It is then required that the forcing amplitude be zero, giving 1 

2 nBq + (,4 — n 2 — q 2 )A = 0. (21) 


4.1 Order 1 


The solution to the first order system of Equation (15) is once again given 
by 


' Xi 


— A x cos(Ar + <p) 

2/i 

= 

kA x sin(Ar + (p) 



A- sin (At + </j) 


1 Of course, the equations could instead be combined into a y /,lf equation, which would 
give a different form of the required relationship. However, under the condition that q 
satisfy the characteristic equation, the two relationships are equivalent. 
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where here the solution is taken as a function of r. 

4.2 Order 2 

Here, in the second order equations, as previously mentioned, x\j n y i/ M and 
z lh appear. Again from the earlier analysis, the linear hub solution may be 
written as 

X\h —A x h cos (A ht + (ph) 

yih = kA xh sm(X h t + <p h ) . (23) 

zih _ A z k si\\{y} x t + Vh) 

Also, the hub frequencies A/> and V} x may be written using a nonlinear asymp- 
totic frequency scaling function where 

&h = 1 + + ld2h 

As already seen, the linear hub frequency is the same as that of the telescope; 
therefore A h — Xuj\ x and v h — i/ujh . For now, also assume that <jJu x — ll>\. 

Typically, a Lindstedt-Poincare analysis is employed for autonomous sys- 
tems. Here, however, the forcing function is an explicit function of time, 
arising from the assumed solution for the hub motion. In the linear huh 
solution of Equation (23), the independent variable t must now be written 
in terms of r. as 

t — t/uj 

T 

1 + LO\ + UJo T ' ‘ ' 

— r(l — uj\ — UJ2 + + • • •■). 

Then, in the hub solution. 

\} x t = Au jfrt 

— A( 1 + U)\ h T ■ * ' )t 

= a(i + + ujoh + * * * )(i — + i — )t 

= A(1 T &2h — -'2 + • ’ • ) T - 

Define the second order contribution to this frequency correction as 

AcJo — UJo h — CJg. 
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where Co represents the frequency correction through order 2. 

Also, write the hub phase angles <ph and ij.'h in terms of the telescope 
angles as 

<ph = 0 + A (f> 

and 

4’h = V + At/;, 

where A</> and At/; represent the offset in phase angles between the hub and 
telescope motion. Substituting into the linear hub solution of Equation (23), 
and retaining the frequency correction through u >2 and ^> 2 / 1 , 

—A x ) j cos(A(l + Acj 2 )t + <f> + A 0) 

= kA xh sin(A(l 4- Aut 2 )r 4- 4> + A</>) . (24) 

A zh sin(^(l - Cu)t + ip + At p) 

This linear hub solution and the first order telescope solution of Equa- 
tion (22) are substituted into the right side of the second order telescope 
equation, Equation (16). After resolving the angles, the three component 
equations are then written as 

x '2 — 2nt/2 — (n 2 + 2A)x2 — — 2uqA(A — nk)A x cos(At + <p ) 

— 3(l 4 ^^BA x A x h 

cos(A(2 T Au-^jr -I- 2 (p 4- A (p) 

— 3(l - B A x Aji, eosfAAuqr + A (p) 

+ | BA z A z u cos((A - i/ h (l - Co))t - At/;) 

— | BA z A z h cos ((A + i/fc( 1 — u>))t + 2t/t 4- Ay;) 

(25a) 

1/2 + 2n.T' 2 + (A — n 2 )t /2 = -2u; 1 A(— k\ + n).4 I sin(Ar + </>) 

3BkA x A x h sin(A(2 -(- Acqj)^" 4- 2</> + A(/>) 

(25b) 
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Co 4 A“ Z 2 — 2lu'i\^ A z sin (At 4 0 ) 

— |S/l x /L/ 1 .sin((A 4- 7/^(1 - a>))r + 04 </; 4 A?/?) 
+ ^BA x A zh sin((A — ^/ t (l — ^’)) r + 0 ~ "0 — A+) 

— ^ B j\ z A x } { sin(A(2 4 Au>>)t 4 i/ ! 4 0 4 Ac>) 

4 sin( A Au^r — 0 T 0 + A0). 

(25c) 

Within Equations (25), each term of the forcing function is treated in 
the manner of the general approach given in Equations (18) and (19). First, 
consider the potentially resonant terms in Equations (25a) and (25b). Using 
the notation of Equation (18). 

41 — — 2ll?iA(A — nk)A x 
B = + n)A x 

q = A 
0 = 0 . 

The condition of Eciuation (21), to avoid resonance terms, gives 

— 2uq A A x [2 ??, ( — k A 4 n) A 4 (,4 - n 2 - A 2 )(A - nk)} = 0. 

In general, this condition can only be satisfied for uj \ — 0. Next, consider the 
resonant-type terms in Equation (25c). Here, 

T - 2uq A 2 /L 
V = A 
7 - 0 - 

The condition to avoid resonance is 

2cjiA 2 -4 = - 0 

and, again, this can only be satisfied for = 0. Accordingly, is now 
taken to be zero. Note that this requirement also gives u = ^ 2 , which will 
be used from this point forward. 

The particular solution to Equations (25) is then built using the method 
of undetermined coefficients, as previously described. The resulting solution 
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takes the form 


x -2 = p 2 \A x A x k cos (A (2 + Aw 2 )r + 2 0 + A0) 

+ PnA x A x h cos(AAw 2 t + Ad)) 

+ P 2 -aA z A zH cos ((A - v h {l - 'jj 2 ))t - A 0) 

+ p2iA z A z f l COS ((A + ///[ ( 1 — Lo’2 ) ) t + 20 + A0) 

'i /2 = sin(A(2 + Au-'oJt + 20 + A0) 

+ 022 A j; A x i, sin(AAw 2 T + A0) 

+ ff23^s^z/iSill((A - !//,.( 1 - 02 2 ))t ~ A0) 

+ o"24' / U^i/i«in((A + - U ; 2 ))r + 20 + A0) 

22 = « 2 i At-^zh sin((A I ///,.(! - it> 2 ))r + 0 + 0 + A0) 

+ K. 2 2/1a^2fcSiu((A - u h {\ - uj 2 ))r + 0 - 0 - A0) 
+ K 2 iA z A xh siu(A(2 + Acj 2 )t + 0 + 0 + A0) 

+ k 2 4j 4 3 A.,./ ( sin(AAta 2 T — 0 + 0 + A 0), 


(26a) 


(26b) 


(26c) 


where the coefficients are 

3B[-2nAJfc(2 + Aca 2 ) + (l + ^)(A 2 (2 + Aw 2 ) 2 + n 2 - A)] 

P 21 ~ A 4 (2 + Aw 2 ) 4 + (A - 2n. 2 )A 2 (2 + Aw 2 ) 2 - (A - n 2 )(n 2 + 2A) 

■3(1 — ^-) B(A 2 Acl) 2 2 + n _ — A) 
p22 ~ A 4 Aw 2 4 + (A - 2?f. 2 )A 2 Aa; 2 2 - (A - n 2 )(n 2 + 2A) 

_ — 2 ^((^ ~ ^(l “ u > 2)) 2 + '»' 2 ~ A ) 

Pz?i (A - Z 2 ft (l - u 2 )) 4 + (A - 2n 2 )(A - i/ h (l - u; 2 )) 2 - (A - n 2 )(n 2 + 2A) 
_ |B((A + v h (l - w 2 )) 2 + n 2 - A) 

/J24 (A + i/ h (l - ta 2 )) 4 + (A - 2;i 2 )(A + u h {l - w 2 )) 2 - (A - n 2 )(n 2 + 2A) 
3 B[ — 2n\(2 + Au-' 2 )(1 + ^-) + k( A 2 (2 + A u 2 2 )" + + 2A)! 

1721 = A 4 (2 + Aw 2 ) 4 + (A - 2?7 2 )A 2 (2 + Aj 2 ) 2 - (A - n 2 )(n 2 + 2 A) 


— 6nAAw 2 (l - ^-)B 

A 4 Ao; 2 ' 1 + (.4 - 2n 2 )A 2 AuV - (A - // 2 )(n 2 + 2/1) 

3n(A - i//,(l - uj 2 ))B 

(A - r/,(l - o> 2 )) 4 + (A - 2n 2 )(A - i/,,(l - a; 2 )) 2 - (A - n 2 )(n 2 + 2A) 

— 3n(A + t//,(l - u,’ 2 ))/j 

(A 4- i//,.(l — ca 2 )) 4 + (A — 2??, 2 ) (A + u h {l — cjo)) 2 — (A — n 2 )(n 2 + 2A) 


71 





3B 


k 21 


^(i 

— ^2)(2A + v/i ( 1 — 

^• 2 )) 




3B 


K 22 


2 ^ 1(1 

-W2)(2A-i/ fc (l- 

o; 2 )) 




3B 


K 23 


2A 2 ( 1 

+ Aa*)(3 + Au/^) 


«24 



3 B 



2 A 2 ( 1 

- Ace 2 ~) 



4.3 Order 3 

The solution procedure at order 3 is much the same as for order 2, with 
much more complicated expressions. The previously determined lower order 
telescope motion is substituted into Equations (17), along with expressions 
for the hub motion. 

For now, the differential equations are examined only for potentially res- 
onant terms — those terms in which A appears as the frequency in the r 
domain. Analysis of these terms will provide relationships between the ac- 
ceptable values of A x and A Z) as well as corresponding expressions for uq?- 
(Further research will allow the development of the actual order-3 solution 
terms.) 

Denoting the Cartesian components of these potentially resonant terms 
by 77 : , : , 7 Z y: and 77-. 

77 ;c — — 2 u> 2 A(A — rik) + 3 C ( ( — 2 )A x i l ~ -f 2L/ 7 ~) 

■h \B[A x h 2 [k(J 2 \ + k°- 22 + 2/921 + 2/9-22) (27a) 

+ A 2 /d(^2i — ^22)] A :v cos(A t 4 - 0) 

7Z y = 2cj’2\(\k — n ) + ~kC A x/t (3A* 2 — 4) + A z ^- 

1 (27b) 

— BA X h 2 ( — kp 2 \ T k P 22 &21 H - ^ 22 ) A x sin (At + <p ) 

77- - [ 2 ^ 2 A 2 + f C((A: 2 - 4)A sh 2 + 3 A zh 2 ) 

— ~B^A Z f^{p 2 3 — pi 1 ) ~ cl.r/ / “(/\23 _ ^21)) + A A z sin(Ar + th). 

(27c) 
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Note that 1Z X and lZ y contain A x as a factor, and 1Z Z contains A z . 

Now examine the conditions to avoid resonance, given by Equation (21), 
and the simpler ^-component condition. It is clear that the coefficients of 
Equations (27) cannot satisfy those conditions with any degree of flexibility 
in the oscillation amplitudes — either or both of /l* and A z would be required 
to be zero. Therefore, it is necessary to pursue an alternative. 

It is possible to increase the flexibility among the orbital parameters by 
introducing additional resonance-inducing terms. One way to do this is to 
require that the hub be in a halo orbit about L 2 . This means that vj x — 
which gives 

1 -^ 2 ) = A(1 + Alc?2 )• (28) 

Making this requirement introduces the following additional resonant terms: 

1Z X ' = 3{C'[cos(Ar — Ad + d + Ad) — cos(Ar + 2d + A d — <p — Ad)] 

+ |5;(fciT23 + ^23 + 2p 2 s) cos(A T — A?/,’ + 0 + Ad) 

+ ( k<j- 2 i + n 24 h 2 ^ 24 ) cos (At + 2 ip + Ad d ~ Ad)]} 

AzA,-.hAzh 

(29a) 

7 Zy = |{C'A:[sin(Ar + 2d + A d - <p — A d) + sin(Ar — Ad + d + Ad)] 

+ 2B[(kp23 — (T 23 ) sin(Ar — Ad + <p + Ad) 

- [kpu + & 2 i) sin (A r + 2 d + Ad - d — Ad)]} 

AzAxh-^zh 

(29b) 

TZj = | {C[(k 2 - 4) sin(Ar + d + Ad - A d) 

+ ( k 2 + 4) sin(Ar - d ~ Ad 4- 2d + Ad)] 

4- 2B[(p 22 - ^"21 ) sin (At 4- d + Ad — Ad) (29c) 

- ( 7^21 + K 22 ) sin (At - d - Ad 4- 2d + Ad)] } 

AxAx-h-Azh- 

It can clearly be seen that the addition of these terms introduces more flex- 
ibility. I 11 particular, note that .1,. appears in 7 Z z ’, while it is not in 77 
similarly, A z is introduced into the x and y direction terms. 

Next, in order to allow these new terms to combine with the original 
resonance terms, it is required that the phase angles match. All in-plane 
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trigonometric arguments are now required to be At + d- out-of-plane argu- 
ments are required to be A t A d- This leads to the following six relationships: 


cos(At A 2'ij) A Ad — d - A0) = A cos( At A 0) (30a) 

cos (At - Aib A d A A<p) = ± cos(At A ©) (30b) 

sin( At A 2d A A d — & — Ad) = A sin(AT 4- 0) (30c) 

sin(AT — Ad A o 4- Ad) — A sin(AT 4- 0) (30d) 

sin(AT 4- 0 A A-y; — A<j>) — ± sin (At -4 VO (30e) 

sin(Ar — y; — Ay; -j- 20 A A<p) = A sin (At -t- d) . ( 3 0 f ) 


From the hub halo motion analysis of [2]. the hub phase angles must 
satisfy the relationship 

0h “ 0h = 2 T 

for arbitrary integer j. With this in mind. Equations (30) load to the re- 
quirement 

d — 0 = (7 A |)tt. (31) 

for arbitrary integer F 

Richardson [2] also indicates that j must be odd (1 or 3) in order to avoid 
enormous hub motion amplitudes which would violate the series expansions 
being employed. Therefore, the same assumption is made here, that j must 
be odd. Applying these relationships to Equations (30) gives 


cos (At A 20 A Ay; — d — Ad) = — ( — 1)* cos(At A d) 

cos(At — Ay; A d A Ad) — ( — 1)^ cos( At A d) 

sin(AT A 20 A A0 — o — Ad) = -(-l) £ sin(AT — <j>) 

sin(AT — Ad A d A Ad) — (— l) f sin(Ar A 0) 

sin (At A 0 — Ad — Ad) = (2 — j) cos(At A d) 
sin(AT — 0 - Ad A 2d A Ad) = —(2 — j) cos(A r A d)- 
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Using these results and augmenting Equations (27) with Equations (29), 


n x 


— 2 o> 2 A(A — nk)A x + 3CA X (k,~'A x h‘ — 2A x f, 2 + A z h 2 ) 

+ t z BA x [(k<7 21 + A:<J22 + 2/921 + 2/922) A-d, 2 + («21 — K22)U-/, 2 ] 

+ ^{—^■) t ^'A i A x hA z h 


+ |(— — fc(J23 + /V23 + 2/923 — ^'0"24 ~ «24 — 2/?24) A Ajft Ah 


cos(A T + cp) 


(32a) 


IZy — 


2uj 2 X(Xk — n)A x + jkCA x (3 k 2 A x h 2 — 4Aft 2 + Ah”) 
+ \BA X A x h ~(—k p 2 I + kp22 — <9"2i + 022 ) 

+ f( — I)' B(kP23 ~ C23 + kp24 + 024) A AhAh 


(32b) 


n z = 


sin (At 4- </>) 

2^2 A~ A + | CA Z (k 2 A x h — 4Ah 2 + 3 Ah 2 ) 

~ f BA, [ ( P23 — P'i-\)A z h — (k<23 — ^ 24) Ah. 2 ] + AA 


- 6(-l ) e CA x A xh A zh 


+ — 1 )* B(p22 + f>2\ + k 22 — ^'21 ) Ac A h -Ah 

(2 - j)(-l/cos(Ar + 0). 


(32c) 


The relationship between Vh and A, given by Equation (28) may now be 
used in the. p, a , and k functions, changing some of their forms to 

-|B(A 2 Au.' 2 2 + n 2 - A) 

Pri = A 4 Aw 2 4 + (A- 2n 2 )A 2 Atn 2 2 — (A — n l ){n 2 + 2A) 

|/A(A 2 (2 + Auj 2 ) 2 + u 2 - A) 

P2t ~ A 4 (2 + Auj 2 ) 4 + (A - 2n 2 )A 2 (2 + Acj 2 ) 2 -(A- n 2 ){n 2 + 2 A) 

_ —AnXAusoB 

° 2:i ~ A 4 Aw 2 4 + (A - 2n 2 )A 2 Aw 2 2 - (/I - n 2 )(n 2 + 2/1) 

— 3?? A (2 -f- z\u.’2)/^ 

CT24 = A 4 (2 + Acj 2 ) 4 + (A - 2n 2 )A 2 (2 + Aoj 2 ) 2 -(A- n 2 )(n 2 + 2 A 
3£ 

K21 ~ 2A 2 (1 + Aa,- 2 )(3 + Aw 2 ) 
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AC 2 2 = 


35 

2 A- (l - Au.> 2 2 ) 


For use in tho equations at third order, the p, a. and k functions need only 
be expanded through linear terms in Au.v When linearized, the functions 
become 

Pi 1 = 8nk\ — (2 + A")(4A' 4- n~ — /l)] 

+ [—^(8 A 2 + .4 - 2n 2 )A 2 [— 8?iAA + (2 + A 2 )(4A 2 + n 2 - /!)] 

+ 4(— nAA + (2 + A 2 )A 2 )] A Ju'2 | 

= p2 l a + PaifcAno 

3(2 - A 2 )/? 

" 22 ~ 2(.e + 2 . 4 ) 

A 

— P22a 

35 

P23 ~ ~2(n 2 + 2.4) 

A 

P 23 1/ 

p 24 = |5{(4A 2 + n 2 -.4) 

- ^A 2 [16A 4 + 8(n 2 - .4) A 2 + (A - n 2 )(3n 2 + .4)]Aw 2 } 

= P f 14a + P'24/)Au.’2 

(72, = |5{[-2nA(2 + A 2 ) + A(4A 2 + n 2 - 2.4)] 

+ [— i(8A 2 + .4 - 2?r 2 )A 2 [— 2nA(2 + A 2 ) + A(4A 2 + rr + 2.4)] 

+ (-n.A(2 + A 2 ) + 4AA 2 )]Aa.’ 2 } 

= CToja + <T 2 u,Aa22 

_ 3?iAA^(2 - A 2 ) 5 
CT22 " (.4 - n 2 )(n 2 + 2.4) 

^ CT22A» Auy’o 

3/iAzXuJ2 13 

aXi = (.4 - n 2 )(n 2 + 2A) 

A A 

a 24 = -fXB{2 - i[48A 4 + 4(.4 - 2n 2 )A 2 - (» 2 - A)( n 2 + 2,l)]Aa> 2 } 

= C 2 4 a + U 2 .l/jAu22 
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B 

«21 = (3 — 4 AW 2 ) 

/\ 

— K.-21o + K2it,Au>2 

ZB 


Kn — - 


K23 


2A 2 

•22c 

B 

6A 2 


a 

— «22a 


(3 - 4Acv 2 ) 


— h -23 o + K23bAuJ2 

ZB 

A 

— «23ai 


where 

d = 16 A 4 + 4(A - 2n 2 )A 2 - [A - n 2 )(n 2 + 24). 

The subscripts 6 and a respectively denote those terms involving and not 
involving Au 2 . 

Keeping this in mind, and recalling that Au 2 is defined as the difference 
between u) 2 h and u>2, Equations (32) may be written in the form 

7Z X = ( Oi I a A;. + O' 1 / ; Co’2 A x + 0:'2«d 3 + COS (A T + 0) 

= (/3l a Ac + ^16^2-4^ + 0 2a Az + ,^26^2-4-) Sill(AT + 0) 

B, z — (yia^4x T T '72a -4 z ~h )(-^ j )( 1) COs(At $), 

where the a, /?, and 7 coefficients are functions of constants and the the hub 
solution parameters u)-2h, Ar/n and A z h- Using a tilde to indicate substitution 
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of CJ2 h f° r Alo’ 2, these functions arc 


a i a — 36' [k 2 A x } } 2 — 24 J / l 2 7 Azh~) 

7 r } B [(A'fJoj 7 Aa 22 7 2/?21 Hr 2p 22 )4 ;r / | 2 7 (^21 “ ^22 ) -4 s / l 2 ] 

O' ib = — 2A(A — nAr) — i^B [(A:<x 2 h, + A '0226 7 2p2ib)A x h 2 + 

n '2a = |( — 1)^ 46 T 7 R{ — k(7-):\ 7 ^23 7 2/9 2 3 — A'YT 2 4 — K'24 “ 2p 2 4 )] 4 x / f ,4-/,. 

0-2 b — “|( — 1)^(~ kaxib + f ''2 36 — ka 2 Ab — ^P'2Ab)Axh^zh 

ft la = |A:C(3A;74 x /7 — 44 x // 2 7 4 2/l 2 ) 7 |B.4. T / f 2 (“A:p2i 7 &p 2 2 ™ &21 + ^ 22 ) 

016 = 2A(AA- — n) + |B4 jr / l “(A;p2it> 7 <7n6 ^ 226 .) 

/^2a = f (— 1/ #(A:/?23 ^23 7 A:p 24 7 ^24 )A T bAzh 

ft‘2b ~ |( — 1)^B ((7 2 36 k p2\b &21b)AxhA-zh 
7la = ( — 1 )^L — 667 7 ~B(p 22 + P21 + ^ 22 ~ ^21 )] Avh^zh 
lib =r 2 ( ^ ) B ( AC 2 1 6 p2 1 b ) -4 a:/i. 4 ; /j. 

72a = j 6 f (fc 2 /l a ./, 2 — 44^7 7 34. /j 2 ) 

7 2 B [(p 2 3 — P24)4 ./j 2 — (k-23 — K, 2 4)At/j 2 ] ~~ A 

72 b — 2A^ 7 |B (p24^4 2 /7 7 «236-4 lX ^ 2 ) - 

( 33 ) 

Recall that the conditions to avoid resonance are as specified by Equa- 
tion (21) and by the requirement that the 2 component of the resonant forcing 
function have zero amplitude. Here, these requirements now take the form 

2n(fti a A x 7 3\bic 2 A x 7 -4a 4- 7 B 2 / 7 2 4 - ) A 

7 (4 — nr — A~)(a j a 4 x 7 7 cv 2a 4, 7 ev 2 / 7 2 4 2 ) = 0 (34) 

and 

7 1 a *4 . 7 ; 7 7 1 / 7 2 4.,. 7 72a 4- 7 72 / 72 4 z = 0. (35) 

Equation (34) may be written in the form 

£i«4 r — c,i7 2-4 x 7 ( 2 a 4, 7 Qib^i^z — 0, (36) 
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where 


Ci a — 2n/?| a A + (.4 — it — A 2 )ai a 

C' 2 a = 2ri/5'2 a A + (.4 — n — A 2 )a' 2a 

T (3/) 

Ci6 = 2n,di&A + (-4 - n - A‘)ai6 

C 2 & = 2n/? 2 &A + (.4 - n - A 2 )a 2 fc. 

Equations (35) and (36) provide two equations in the three variables A K , 
Aj, and wo. However, a third relationship may be introduced. Consider the 
order- 1 solution of Equation (22). Using the frequency and phase relation- 
ships of this section gives 


Si 


- A x cos(Ar + < j >) 

V \ 

= 

k /\ L . sin(Ar — < f >) 

_ z \ 


_ (2 - j)(-l)C4 £ cos(Ar -4- < j >) _ 


with j either 1 or 3, and £ either 0 or 1. This part of the solution corresponds 
to the case where the aperture plane maintains a roll angle £ about the y- 
axis. The range of this roll angle is related to the integers j and £ in the 
fashion shown in Table 2. 

Table 2: Relationship of Aperture Plane Roll Angle to i and j 


j 

£ 

£ 

1 

0 

(0, 7t/2) 

1 

1 

(-tt/2,0) 

3 

0 

(-tt/2.0) 

3 

1 

(0, tt/2). 


Therefore, the ratio (2—j)(—l) t A z )/A x may betaken as an approximation 
to tan £. Let q represent A z /A x . Assuming A x to be nonzero, Equations (35) 
and (36) may be written in terms of q as 

7la ~ 716^2 + 72a 7 + 726^2'/ = 0. (38) 

and 

Ci« + Cua J, 2 + C'ia q + Qzh^iq = 0 (39) 
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These two equations may he solved simultaneously for <jj 2 and //. The 
combination of the equations leads to a quadratic, with two sets of solutions 
for lj-j and rj. Depending on the hub motion amplitudes, the solutions for 
r/ correspond to approximate aperture plane roll angles from —n /2 to n/ 2 . 
The solution for jjo is then 

u* = (40) 


where 


A = jibQib — 72&C16 

& — 7 1 a Qlb T 7 ih^2a f'2a C 1 b /2/Ala 

C 7laC2a 72a Cl a* 


and the corresponding q is 


(la + (l/A ? 2 
Qla + 


r I a 


7 1/^2 


■f 2a 


yob ^2 


(41) 


In turn, each value of // is associated with an infinite set of (A J: .A Z ) pairs, 
each of which is, in turn, associated with a set of initial conditions for the 
relative telescope motion. 

A rule of thumb may be developed in order to bound the selection ol A x 
and A z . Examining the solution to the linear telescope equations, given by 
Equation ( 22 ), it can be seen that the approximate maximum distance of a 
telescope from the hub is 


dm ax = max( \/A x 2 + A z \ kA x ) 
= A x max( \J 1 — rf, k). 


(42) 

(43) 


Therefore, if the maximum distance is required to be no greater than a dis- 
tance d max , A v may be selected such that 


A x < 


(in 


max( y/l + rf. k) 


and then A 2 = 77 A r . 
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4.4 Solution Summary 

The solution is now constructed by adding Equations (22) and (26), and using 
the frequency and phase angle relationships of Section 4.3. This combined 
solution through order 2 is then 

x = —A x cos (At 4 d>) 

4 (p 2 \A x A x h — p 2 iA z A Z h{—\) t ') eos(A(2 + lS.lo 2 )t + 6h 4 <J») (44a) 

+ {p 7 nA x A xh 4 • P 2 3 A z A 2h (-lY) cos(AAw 2 t + 6 h - <p) 
y = kA x sin(Ar 4 #>) 

+ (cr 2 i A x A xh - (T 24 A z A zh (-lY) sin(A(2 A lj 2 )t + 4>h + <t>) (44b) 

+ [cr 2 2 A x A xh - <J 2 zA 2 A th (-lY) sin(AAu> 2 T + o h - 6) 
z= [(-l) f /Ucos(Ar + o) 

+ (xviAxAzh 4 K 2 :iA z A x h(— 1)^) cos(A(2 4 Ac o 2 )t 4 (frh 4 <p) 

- ( K<aA x A zh 4 K 2i A z A xh (-l) e ) cos(AA u 2 t + <Ph - (p)} (2 - j), 

(44c) 

where r = (l + oj 2 t). The coefficients are as given following Equations (26), 
with the modifications which follow Equations (32). 


5 Simulation 

5.1 Implementation Procedure 

Various representations of the telescope equations of motion were coded in 
MATLAB to examine the differences among the solutions. Some of the results 
are presented in this section. The motions of the hub and only one telescope 
wore simulated because of programming and interpretation simplicity at this 
point. In the future, a. visualization could be developed based on selection of 
initial conditions and results from the simulations. 

The selection of the initial state is not arbitrary because the resulting 
uncontrolled orbit path is likely to be unsatisfactory. The solutions of these 
models show that the hub’s trajectory diverges from a closed, periodic orbit 
in typically 150 to 200 days. Exploration of initial conditions has been a time 
consuming portion of this research. 
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Additionally, the uncontrolled telescope motion will diverge from the 
nearby reference path about the hub. 

Difficulties in finding initial conditions that produced useful trajectories 
motivated short duration time frames. Assume that as the aperture plane 
orbits the L 2 point, a typical observation will last only a few days, before 
reorienting the plane to observe some other target. By simulating runs of 20 
days duration, small diverging motions from the nominal halo motion were 
maintained. 

The baseline solution is the numerical integration of the full equations 
of the huh subtracted from that of the telescope, as shown in Equation (1). 
In essence, there are two separate vehicles orbiting the L 2 point with very 
similar initial conditions. The slight differences in initial conditions are due 
to the placement of the telescope with respect to the hub. 

One cannot arbitrarily specify initial hub and telescope positions for either 
the full, or the truncated nonlinear model, or the analytical model, and obtain 
closed or nearly closed orbits about the L 2 point. 

To begin the process, implement the lengthy algorithms shown in Sec- 
tion 4.3: 

• begin with a selection of the amplitude of the motion along the 
^-axis of the hub with respect to L 2 

• use a consistent set of physical quantities taken from Dunham and 
Muhonen [3] for the constants used in computing A , B : C of Equa- 
tions (8) and for A and k 

• compute the o:. J. 7 terms of Equations (33) 

• compute the (," terms from Equations (37) 

• solve for both of the : x 2 solutions from Equation (40) 

• maintain two solutions for 7/ (Equation (41)) based on two terms 

• then A z — 77 A x yields two choices 

Table 3 shows three combinations of amplitude selection, selections of j. 
and the results of calculations based on these selections. The application 
of Equation (G-l) from Richardson ;2] constrains the A x f, corresponding to 
the selected A z } x . (This equation is plotted in Richardson's Figure' H-l.) 
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Table 3: Summary of Inputs and Resulting Values Used to Select Orbit Size 
and Orientation 


a x /i 
(km) 

(kin) 

j 

U'2 

■jJ 2 + 

9 

n + 

r 

(deg) 

V 

(deg) 

A-~ 

(m) 

A s + 
(km) 

‘280,007 

050.150 

1 

-0.20211 

-0.20185 

1.6115 

1 .6147 

-58.18 

-58.23 

80.58 

0.08073 

280,667 

550.459 

3 

-0.29211 

-0.29185 

1.6115 

1.6147 

58.18 

58.23 

80.58 

0.08073 

227,219 

250.000 

1 

-0.26069 

-0.11043 | 

0.45009 

3.9851 

-24.23 

-75.92 

22.50 

0.1993 

227,210 

250.000 

3 i 

-0.26060 

-0.11043 > 

0.45000 

3.0851 

24.23 

75.02 

22.50 

0.1993 

211,126 

1,000 

1 

-0.23166 

-0.07373 

0.00178 

905.98 

-0.1020 

-89.94 

0.08902 

45.30 

211,126 

1.000 

3 

-0.23166 

-0.07373 

0.00178 

905.98 

0.1020 

89.94 

: 0.08902 

45.30 


The amplitude of the hub motion along the 2 -axis was first chosen and then 
the amplitude of the hub motion along the x-axis is constrained by this 
relationship, which assures halo orbit solution via correct frequency selection. 

The choice of .4./, of 550,459 km is the largest z- value that can be selected 
and yield real roots for ui 2 - An intermediate value of 250,000 km and a very 
small value of 1,000 km were also selected for presentation. The appendix 
shows the results of many calculations of this example. 

There are only two acceptable choices for j, which controls the sign on 
the z-position equation. The selection of j = 1 corresponds to a left roll 
orientation of the aperture plane; the selection of j = 3 corresponds to a 
right roll; This table shows the results of calculations based on fixing 6 = 1 
and A x = 50m. (Results for 6 = 0 are not shown because the q values would 
be negative in Table 3 and lead to negative values of A-.) 

Because the equation for u 2 admits two roots, these two values are carried 
forward to see their effects upon the calculation for q. The superscripts — 
and + refer to the sign on the square loot in each solution. Two values are 
calculated for q and, in turn, two values of the arctangent of q to yield the 
approximate aperture plane roll angle, £. 

The resulting magnitudes of the q terms, and their corresponding angle 
values should be highlighted. As stated above, the large value of A z & is 
associated with the case in which the two roots for u >2 approach a single 
double root. Accordingly, the resulting values of q, f, and ,4, are also nearly 
identical. 

For the intermediate value of A-/,., the u >2 values are different. The roll 
angles, of the two solutions, are quite different. Finally, at the smallest 
value of A, fl . the u >2 solutions are quite different. The roll angles of the two 
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solutions arc at extreme values based on rj ~ , the halo plane is nearly in 
the ecliptic (x-y) plane and based on //+, the halo plane is nearly in the y-z 
plane. 

Because the user can select either A z ^ or the desired roll orientation, 
the user does not have control over the combination of halo plane size and 
orientation. When considering the roll orientation of the aperture plane, 
from — tt/ 2 to tt/ 2, the user should realize that three months later, the orbit 
will have moved about the Sun such that the initial roll angle appears as a 
pitch rotation. This offers some orientation flexibility. 

Again, the user does not have total control over hub position and. of 
course, the aperture plane location with respect to L 2 . 

Continue with the implementation of these algorithms: 

• compute the p. a, k. terms on pages 76 - 77 

• compute the analytical expressions of Equations (44) 

Equations (44) are the analytical expressions for the motion of the tele- 
scope with respect to the hub. By taking the derivatives of these equations, 
one can compute the full state as a function of time. These equations may 
be used by trajectory control designers. 

Also at t = 0, Equations (44) may be used to determine the initial state 
for use by two comparison solutions requiring numerical integration: 

• Solution of Equation (1) - full nonlinear equations of motion of the 
telescope with respect to the hub 

• Solution of Equation (8) - truncated, to second order in the distance of 
the hub from /> 2 , nonlinear equations of motion of the telescope with 
respect to the hub 

5.2 Simulation to Compare Results of Models 

Here is an example that compares the results of the three solution methods 
for calculating the motion of the telescope with respect to the hub. 

: ‘FulP refers to the numerical integration of the separate equations of 
motion for the telescope and the hub, each spacecraft with respect to L 2 , and 
then subtracting the hub motion from that of the telescope. This represents 
the implementation of Equation (1). 


Table 4: Initial Conditions of Telescope With Respect to Hub Resulting From 
Analytical Solution 


*(0) 
y( o) 
*(0) 

— 64.7796 m 
0.0 m 

26.4452 m 

±(0) 

m 

~(0) 

0.0 m/day 

4.4205 m/day 
0.0 m/day 


"Truncated” refers to the numerical integration of the truncated equa- 
tions of motion for the telescope with respect to the hub. This represents 
the implementation of Equation (8) with the terms that have a subscript h 
referring to the linear expressions for the motion of the hub. The linear hub 
solution is given by Equation (23), which was further simplified to give a 
halo orbit setting both frequencies to the same value and setting both phase 
angles to zero. 

The initial conditions for both the full and truncated numerical integra- 
tions are obtained from the calculation, at t = 0, of the analytical solution. 

"Analytical” refers to the second order solution, which is an explicit func- 
tion of time, given by Equations (44). 

From the intermediate case of Table 3, select the out-of-plane hub am- 
plitude from the Lo point to be 250,000 km (A z h). Using Equation (G-l) of 
Richardson [2], the resulting value of A xh ■ is 227,219.42 km. Select the row 
with j = 1 or 3, corresponding to the desired aperture roll direction. The 
roll direction is implemented in the solution via Equation (44c). The roll 
angle direction is also indicated in Table 2. Remembering that Table 3 was 
computed using ( = 1, choose j — l for a negative roll or choose j — 3 for a 
positive roll. In this simulation, j = 3. 

Select the initial amplitude between the hub and the telescope to be 
A x = 50 m, along the x-axis. The A z ~ choice initially gives the telescope a 
linear amplitude of 22.5 in away from the hub along the z- axis. 

The analytical solution was first evaluated at t = 0 to obtain the initial 
conditions (of the telescope relative to the hub, Table 4) for the full and 
truncated solutions. 

The initial conditions for the hub are listed in Table 5 with the equations 
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Table 5: Initial Conditions of Hub With Respect to Sun-Earth L 2 Point 


ffc(0) 

227.219.42 km 

yi,( o) 

0.00 km 

•3,(0) 

250,000.00 km 

0 

0.00 rad 

<Ph 

0.00 rad 

ifc(0) 

( A/A:)y,(0) km/day 

?),(0) 

— kXxh (0) km/day 

2,(0) 

-(2 - j)\A zh xin((p h ) km/day 


for x 7 y , and z reproduced here for convenience. These were given earlier in 
Section 3.3. 

The complete initial conditions for the telescope with respect to the L 2 
point are given below: 

*t(0) = .^(0) + x{ 0) 

yd 0) - Vk( ;o) + y(o) 

*t(0) — 2h{ 0) + z( 0) 
i*(0) = ;tfc(0) + x{0) 
yd 0) - Vh{0) + 2/(0) 

Ct(0) = 3?/t(0) + i(0) 

The following plots show the small differences among the three solution 
types over a 20-day duration. The plots show the telescope motion with re- 
spect to the hub. The origin of the coordinate system represents the position 
of the hub. Scales are expanded for the <r and c positions. 

Returning to the discussion that began in the introduction to Section 5.1, 
we note that over 5 days each solution yields almost the same result. The 
analytical solution is quite accurate for trajectory control analysis. We as- 
sume that as the aperture plane orbits the L 2 point, a typical observation 
will last only a few days, before reorienting the plane to observe some other 
target. By simulating runs of 20 days duration, we can see small diverging* 
motions from the reference numerical solution. 

In Figure 3, observe that both the full and truncated solutions are essen- 
tially the same. Considering the scale of the axis labeled ;t x position”, the 
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y position (m) ^ x position (m) 


-60 



Full Truncated Analytical Trunc C=0 


•o 3: Compare 3 Solutions for Telescope Motion Along the .r-axis. 



Figure 4: Compare 3 Solutions for Telescope Motion Along the y-axis. 
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Figure 5: Compare 3 Solutions for Telescope Motion Along the 2 - axis. 

analytical expression is close to the other solutions over 5 days. The curve 
labeled "Trune C=(T shows that when the lowest order term of the truncated 
solution is set equal to zero, the effects are not significant in the ^-direction. 

In Figure 4, observe that all solutions are essentially the same for the 
scale shown. 

In Figure 5, we see that both the full and truncated solutions are essen- 
tially the same. The analytical solution stays close to these solutions over 
5 days. Remember the truncated solution has all of terms A : B, and C in- 
cluded. The highlight is the curve labeled u Trunc C=0” showing nearly the 
same result as the analytical solution. Now the importance of the Ci C'' term 
is evident. It indicates that a more complete analytic solution would include 
higher order terms. 


6 Summary 

This report details the preliminary work describing the formation flying be- 
tween spacecraft near the Sun-Earth L 2 libration point, beginning with the 
circular restricted three- body problem for the hub motion about L 2 . 
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The halo orbit was specified as a desirable aperture plane orbit, instead 
of a Lissajous orbit, because it is periodic. Additionally, over a one or two 
day observation period, the spacecraft telescopes of the aperture plane will 
minimally separate due to natural motions. 

The analytical solution, Equations (44), for the motion of a typical tele- 
scope relative to the hub is presented, including terms which are linear in the 
hub motion. We anticipate feedback from our sponsor concerning the utility 
of these equations for use in orbit control system design. 

We developed a full non-linear solution of hub motion, with respect to 
L 2 , subtracted from telescope motion, with respect to Z/ 2 , as a reference 
trajectory. 

We developed a truncated non-linear solution to the telescope motion 
with respect to the hub. Both the non-linear model and the truncated model 
are solved by numerical integration. 

All three models produce trajectories that diverge from a closed path over 
the course of one orbit (roughly 200 days). We compared the three models 
over a shorter duration of 20 days expecting that a science observation would 
be on the order of a few days, before reorienting the aperture plane. 

The analytical solution algorithm is quite lengthy. One example was 
presented and used to compare the accuracy of the solution compared to the 
reference trajectory. Over 5 days, the analytical solution is very close to the 
reference trajectory. 
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Appendices 

A Constant Parameters 


Table 6: Physical Constants [3] 


Pi 

132,712.440,017.987 km 3 /s 2 

P2 

403,503.236 km 3 /s 2 

n 

0.199106385 xlO -0 rad/s 


Table 7: Derived Constants 



151,105,099.094445 km 

£>i 

454.84086785372 km 

D 2 

149,597,415.850132 km 


Table 8: Constants Specific to Sun-Earth L 2 Point [2] 



-14.82882087 

1/(DU 2 -TU 2 ) 

^2 

1.673691389 

1/(DU 2 -TU’ 2 ) 

5 2 

-0.7444513767 

1/(DU-TU 2 ) 

•52 

0.1250471558 

1/(DU-TU 2 ) 

DU 

1.507683382 x 10° 

km 

TU 

58.13235527 

days 


Richardson’s algorithm ([2], page 2-31) for computing A x h given A z h and 
for computing u> 2 h are given here: 

+ A — 0 
^'2 h = *1 A]. h + 


91 






Those arc the modifications used to dirnensionalize Richardson’s algo- 
rithms: 

A xh = DUy/(-£2(A zh /DU) 2 - ATIPWi 
U>u, = sMzh/DU) 2 + s 2 (A zh /DU) 2 


Where in Table 8. the constants ;i DU” refers to distance unit and "TU” 
refers to time unit. 
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Tabic 9: Computed Values and Coefficients 


A 

B 

C 

1.16605228517927 x 10' ;j 1/dav 2 
5.84853993419721 x 10" 10 l/(km-day 2 ) 
3.86667873919725 x 10" le l/(km 2 -da.y 2 ) 


3.53850956958284 x 10" ? rad/day 
3.18712225987377 

1.16605228517927 x 10“ 3 rad/day 2 
8.60527122236636 x 10~ 5 1/dav 2 

d 

P2la 

P2V, 

P22a 

P22b 

P2'ia 

P2'ib 

P24a 

7*246 

2.56733055729898 x 10“ 5 1/dav 4 
1.18887100881492 x 10 -6 1/km - 
-6.40826960245161 x 10~ 7 1/km 
-2.72318375685819 x 10" 6 1/km 
0 1/km 

-3.33815613931628 x 10“ 7 1/km 
0 1/km 

1.41409729921967 x 10" 7 1/km 
-1.21027861225803 x 10- 7 1/km 

&2la 
0216 
&22a 
&22b 
0 r 23a 
17 24b 
024a 
&24b 

6.51772783060505 x 10" 7 1/km 
-7.61518233174130 x 10~ 7 1/km 
0 1 /km 

—3.81021051605070 x 10 -6 1/km 
0 1 /km 

4.67066447286556 x 10~ 7 1/km 
-8.32024709778584 x 10~ 8 1/km 
1.30305530510727 x IQ" 7 1/km 

K 2l a 

tf'21 b 

*'22 a 

K 22 b 
k 23 a 

«23 b 

^24 a 
^•246 

2.33548302511691 x 10“ 7 1/km 
-3.11397736682255 x 10" 7 1/km 
7.00644907535074 x 10~ 7 1/km 
0 1 /km 

2.33548302511691 x 10“ 7 1/km 
-3.11397736682255 x 10“ 7 1/km 
7.00644907535074 x 10~ 7 1/km 
0 1 /km 
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Part 3: Modelling the Perturbations 
— Elliptical Earth Orbit, Lunar 
Gravity, Solar Radiation Pressure, 

Thrusters 


97 



Contents 


1 Introduction and Problem Definition 105 

2 Equations of Motion 109 

2.1 Circular Restricted Three-Body Problem 109 

2.2 Elliptical Restricted Three-Body Problem Ill 

2.2.1 Elliptical Problem Libration Point Analog Ill 

2.2.2 Relative Motion Equation Derivation 112 

2.3 Lunar Gravitational Effects 116 

2.4 Solar Radiation Pressure Effects 117 

2.5 Spacecraft Thruster Effects . . . . 118 

2.6 Summary 121 

3 Expanded Equations of Motion 126 

3.1 Circular Restricted Three-Body Problem 126 

3.2 Elliptical Restricted Three-Bodv Problem 126 

3.3 Lunar Gravitational Effects • 129 

3.4 Summary 134 

4 Modeling Uncertainties 136 

4.1 Elliptical Restricted Three-Bodv Problem 136 

4.2 Lunar Gravitational Effects 137 

5 Sensitivity Summary 138 

6 Summary and Conclusion 139 

Appendices 143 

A Sensitivity of Telescope Motion to Errors in Hub Position 143 

A.l Derivation of Variational Equations 143 

A. 2 Linear Equations 147 

A. 3 Solution Development 148 

A. 3.1 Non-secular Solution . . 149 

A. 3. 2 Secular Solution 150 

A. 3. 3 State Transition Matrix 153 

A. 4 Out-of-Plane Solution 153 

A. 4.1 Secular 154 


99 



A. 4.2 Non-Secular 156 

A. 5 In- Plane Particular Solution 157 

A. 5.1 Secular 158 

A. 5.2 Non-Secular 161 

A. 6 Summary 163 

B Relative Motion Near Earth-Moon L 2 Libration Point 164 

C Calculation of Luminosity Reduction Due to Partial Eclipse 165 

D Constant Parameters 171 

References 172 


100 



List of Tables 


1 Acceleration due to Solar Radiation Pressure vs. Surface Area 122 


2 Example Initial Conditions 124 

3 Along-Ellipsoid Effect of Sun-Earth Perturbation on Solution 

(90 days) 127 

4 Along-Ellipsoid Effect of Lunar Perturbation on Solution (90 

days) 133 

5 Along-Ellipsoid Effect of Earth-Moon Perturbation on Solu- 
tion (15 days) 164 

6 Physical Constants 171 

7 Computed Values and Coefficients 171 


101 


List of Figures 

1 Aperture Formation of MAXIM Pathfinder, Concept of June 

2002 ' 106 

2 Coordinate Axis Definition 109 

3 Coordinate Axis Definition (including Moon) 117 

4 Vehicle Thrust in Body-Fixed Frame 119 

5 Frame Rotation from Body to Rotating Frames 121 

6 Effects of Perturbations on Relative Distance — Full Equations 125 

7 Effects of Perturbations on Relative Distance — Expanded 

Equations 135 

8 Relative Motion Error Caused by 1.7 km Initial Hub Position 

Error 139 

9 Relative Motion Error Caused by 17.3 km Initial Hub Position 

Error 140 

10 dz/dxf, and dz/dz), with resonance 156 

11 dz/dzh without resonance 157 

12 dx/dyh and dy/dyn with resonance 160 

13 dx/d:Ch and dy/dxf,. with resonance 161 

14 dx/dzh and dy/dz h with resonance 162 

15 dx/dz h and dy/dz h without resonance 162 

16 Solar Radiation Shadowing Geometry 165 

17 Earth Obscuration and Solar Disk Geometiy 166 


103 



1 Introduction and Problem Definition 


Some continuing long-term goals of our sponsor are to develop high-fidelity 
equations of motion representing the formation flying of a spacecraft con- 
stellation near the Sun-Earth L 2 point and equations describing the relative 
motion between these constellation members. 

This is our third report, which continues from Part 1 [1] and Part 2 [2]. 
In this report, we further the work of Part 2 by developing the elliptical 
restricted three-body problem from the previous work with the circular re- 
stricted three-body problem. Additionally, this new work includes the force 
perturbations of lunar gravitation, solar radiation, and spacecraft thrusters. 
The effects of the elliptical orbit of the Earth-Moon system about the sun 
and the force perturbations are incorporated as additive perturbations to the 
circular restricted three-body problem. 

This work builds from the previously developed circular restricted three- 
body problem formulation. The familiarity of the formulations gained from 
Part 2 is maintained here. We present the development of the full nonlinear 
baseline model, which includes these perturbations: 

• begin with circular restricted throe-body problem 

• then adding terms describing the elliptical orbit of Earth around the 
sun 

• then adding terms describing the moon’s motion about Earth 

• then adding terms incorporating the force due to solar radiation 

• finally adding terms for spacecraft thrust 

We identify the magnitude of the contribution of each perturbation to 
the solution, in order to help determine when it should be included in the 
computations. We identify substantial modeling uncertainties. In addition 
to identifying the magnitudes of the effects of the included perturbations, we 
identify the magnitudes of what is dropped in the truncation process. 

■ The development continues with the expanded form of the full nonlinear 
baseline to sufficiently high order such as to capture the relevant contribu- 
tions to the model. The presentation of the perturbations is organized in 
similar fashion to the baseline presentation so the reader can easily compare 
and contrast the models. 
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This research is loosely motivated by formation flying concepts for the 
MicroArcsecond X-ray Imaging Mission (MAXIM) Pathfinder. A concept 
configuration of the MAXIM mission is depicted and described below. 



Figure 1: Aperture Formation of MAXIM Pathfinder. Concept of June 2002. 


Seven spacecraft shown in Figure 1 are in the same plane forming a flat 
aperture. The optical hub is in the “center 1 of six free-flying spacecraft 
(telescopes). Dashed lines do not indicate a material connection, but rather 
indicate random radial directions listed r\ through . The scalar length of 
r ranges from 100 to 500 meters. Additionally, the six free-flying spacecraft 
are loosely spaced 60° from one another. 

There is also a detector spacecraft located approximately 20,000 km and 
90° out. of the plane formed by the flat aperture. The entire system is in a 
Sun-Earth L 2 orbit, which has yet to be designed by the mission team. The 
analysis of the detector’s orbit relative to the aperture plane is beyond the 
scope of this work. 
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In Section 2.1, we present the circular restricted three-body problem’s 
equations of motion, the coordinate frame that describes the motion, and 
define terms used in the formulation. 

In each section below, as appropriate, we explain the effects of those terms 
retained and dropped from the derivation in the course of the truncation 
process. 

• In Section 2.2, we present the elliptical restricted three-body problem, 
which is used as the baseline for the subsequent development., 

• In Section 2.3, we incorporate the effects of the lunar gravity. 

• In Section 2.4, we incorporate the effects of the solar radiation pressure. 

• In Section 2.5, we incorporate the terms used to include spacecraft 
thrusters. 


• In Section 2.6, we consolidate the equations into one location. There 
is further discussion to compare the relative magnitude of the terms so 
that they may be turned on or off. 

• In Section 3.1, we present the expanded form of the circular full non- 
linear baseline. 

• In Section 3.2, we present the derivation of the expanded form of the 
elliptical full nonlinear baseline. 

• In Section 3.3, wo incorporate the effects of the lunar gravity. 

• In Section 3.4, we consolidate the equations into one location. There 
is further discussion to compare the relative magnitude of the terms so 
that they may be turned on or off. 

Because there is no expansion of the forces from solar radiation pressure 
and spacecraft thrust, the contributions presented in Sections 2.4 and 2.5 are 
merely restated in Section 3.4. 

One important measurement of the hub-telescope concept is knowledge 
of the distance between the spacecraft to within millimeters or even smaller. 
However, while many of the equations of motion presented in Sections 2 and 3 
have terms that require numerical values which are listed in readily available 
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publications, some values of the needed physical parameters are given to a 
number of significant figures with uncertainties in their knowledge or mea- 
surement. A certain value may be defined one way in one text and a different 
way in another text. This will produce different results in the numerical cal- 
culations. Examples of these uncertainties applicable to the equations in this 
report are shown in Section 4. Some discussion is provided indicating the 
impact to the simulation results due to differences in the calculations. 

Section 5 summarizes the results of Appendix A. which contains a detailed 
discussion of the sensitivity of the relative motion to small errors in the 
assumed position of the hub. 

Appendix B fulfills a request by the sponsor to minimally address relative 
motion at the L z point of the Earth-Moon system. 

Appendix C explains the eclipse geometry and corresponding shadow 
model used in the application of solar radiation pressure. 
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2 Equations of Motion 

2.1 Circular Restricted Three-Body Problem 

In Parts 1 and 2 of this research [1. 2], the general second order differential 
equations of motion were constructed for an object near the Sun-Earth Li 
libration point, using the force model of the classical circular restricted three- 
body problem. In this model, Earth is treated as being in a circular orbit 
about the sun, the spacecraft mass is considered to be negligible as compared 
to the two primaries, and only point-mass gravitational forces are considered. 


▲ 

A 

y 



Figure 2: Coordinate Axis Definition 

For this system, depicted in Figure 2, the differential equations of motion 
for an object (object i) near the Sun-Earth Li are given by 

f ( = _ (A_ + j rj . ( '' - ■(? ; + g if + * + a,,*, 

VPH 3 P 2 , J \ PlC P 21 )' 
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where 


r, ; — vector from L 2 to object i 
H i = solar Keplerian constant 

//. 2 — terrestrial Keplerian constant (Earth + Moon) 

PH - distance from Sun to object i 
P'2i = distance from Earth-Moon barycenter to object i 
x e — distance from system barycenter to Lo 
D\ - - distance from system barycenter to Sun 
Dy — distance from system barycenter to Earth-Moon barycenter 
x = unit vector parallel to Sun-Earth line of syzygy, 
pointing in Sun-to-Earth direction 
a = terrestrial mean motion about Sun (assumed constant). 

The coordinate frame of Figure 2 is a rotating reference frame with origin 
O at the system barycenter. The x-axis points along the Sun-Earth line of 
syzygy and the c-axis is parallel to the system angular momentum: the y-axis 
completes the dextral coordinate system. 

Let rf r and r* denote, respectively, the vector from L 2 to the hub and to 
a telescope. Therefore, if r is the vector from the hub to the telescope, the 
differential equation of motion for the telescope relative to the hub is 


r = r t - r h 


P\ ^2 

Pit 3 P 2 t 3 


r t 


Pi 


+ 


P'i 


Pih 3 Pih 3 


r h 


/q(x e + Pi) fi 2 (x e - P 2 ) 
Pit 3 pit 3 

fi x {x e + D\) /q(x, - D->) 


-/q 


r/ 




Pit 3 Pih 3 


- 1*2 


Pih . 3 

r t r h 


PvP 


X 


P 2 3 P 2 I 3 


( 1 ) 


T D\) ( 3 3 ) ^ D‘ 2 , ) ( 3 3 

K Plt Pih ) \Plt P2h 


X. 


where now, in general, the subscripts h and t refer to the hub and telescope. 
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2.2 Elliptical Restricted Three-Body Problem 

The extension to Equation (1) is now developed for the case of the elliptical 
restricted three-bodv problem. First, it is necessary to locate the point which 
is analogous to the libration point L 2 in the circular restricted problem. As 
would be expected, such a point exists; as the Sun-Earth distance varies, the 
location of the point oscillates along the x-axis. 

2.2.1 Elliptical Problem Libration Point Analog 

Say that there is an elliptical analog to the L% point and that its position 
relative to the (assumed inertial) system barveenter is given by 

R e = x e x + y e y -+- z e z. 

Letting / refer to the true anomaly of Earth’s orbit about the sun, the 
coordinate system has angular velocity 

= /z, 

which now is considered to be varying throughout the year. Differentiating 

Re, 

ie “ he 

Ve + fXe i 

where the column vector notation is used to indicate the xyz vector compo- 
nents. 

The Newtonian gravitational force per unit mass acting on an object at 
this point is given by 

Pl{Xe + Di) fl 2{%e “ D 2 ) 

fMVc _ l>l2Vc 

fie 3 »’2e 3 
_Ml£e _ 

J”le 3 f'2e 3 

where 

rie = {x e + £>i)x + y e y + z e z and r 2e = {x e - D 2 )x + y e y + z e z. 





%e file 2 f y e j X e 

ij e + + 2 fx e - f 2 y e 

i e 
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and D\ and Do now refer to tho time-varying locations of Sun and Earth 
along the line of syzvgy. 

Clearly, as in the circular restricted case, the 2: equation is decoupled, 
and admits a solution z v -- 0. If, as anticipated, the desired solution involves 
y e — 0, the x and y components of the force- acceleration equation become: 


x : x e - f 2 x e + 


hi 


+ 


h 2 


(:r f . + D\ ) 2 (x c — D 2 


- 0 


y : /. Ct +2/x e =i|(x t 2 /)=0 


From the y-coinponent equation, x e 2 f is therefore constant. However, con- 
servation of the Sun-Earth two-body angular momentum per unit mass h, 
gives 

D\f - K 

where D is the varying Sun-Earth distance D\ + Do- Of course, unlike the 
case of the circular restricted problem, this distance varies throughout the 
year. Substitution for / in the y ~ component equation gives 



constant. 


For later convenience, define this constant in terms of the constant 7 , such 
that 



Taking the positive root, this gives the definition of 7 as 


A Ze ~ D 2 

7 “ D 


In the case of the Sun-Earth system, 7 is approximately 0.01007824: x e 
varies between 1.486 x 10 8 and 1.536 x 10 s kin throughout the course of the 
year. 


2 .2.2 Relative Motion Equation Derivation 

Let R ? denote the position vector from the system barycenter (point O) to 
spacecraft ?n L . Referring again to Figure 2. 

Ri = r c x f iy. 
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As above, the coordinate system has angular velocity 


- fz. 

Differentiating R,, 

R t = x e x + fx e y + r, 

Rt = (Ze - f 2 X e )x + (2 fx e + fx e ) y + l'i. 

The distance x ( , and its time derivatives may be written in terms of the 
varying Sun-Earth distance D and its derivatives, using the following defini- 
tions of the constants 7 and p, and the associated relationships: 

D x = pD 

D 2 = (1 - p)D, 


= (7 + 1 - p)D 
x E = (7 + 1 - p)D 
x ( . = (7 + 1 - p)D. 

This gives the acceleration vector as 

R,; = (7 + 1 - p)(D - f 2 D)x + (7 + 1 - p)(2fD + fD ) y + f.,. (2) 

Similarly, two-body relationships may be used to write D and D in terms 
of / and /. Using the two-body equations of motion for the sun (could use 


7 + 1 — 


A x e - Do 
D 

X tr + D\ 


D 


A 

P = — 
P 

l-P-*- 

P 


where 


Then, 


P = Pi + P2- 


X e — 7 D + Do 
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Earth): 


Ri = — £>ix 

= —pDx 

Wc = /Z 

R, = —pDx. - pj'Dy 

R, = {-pD + pf z D)x + (-2 pfD - pfD) y 
= P(J" 2 D - £>)x - pU'D + 2fD)y 
— -j^x (two-body force) 


This gives the component equations 

D 2 


x : p(f 2 D -D) = ^ 


y : -p(fD + 2fD) = -^±(D 2 f) = 0, 


which vield 


D - f 2 D = 


/<2 

pD 2 


2/1) + fl) = 0. 

Substituting into Equation (2), the acceleration becomes 

= -(7 + 1 - + fj 

■(7 + 1)” + 7TT I x — E 


-(7+1)0 


pD 2 

H ( /I o 

D 2 + D 2 

P-> 

D 2 


x + r,-. 


Introducing the two-body forces, 

R,: = 


/ ■ i i \ /'i „ /'2 

(7 + f)^- f D 2j 


x + E = -- ^pii - -^P2r (3) 

P\i P2i 


Note that the explicit appearance of the time derivatives of / has now been 
removed. 
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Next, the vectors p 1( and p 2i may be written as 


Pit — ( X e + Di)x + Ti 
= (7 + l)Dx + r,; 

Substituting into Equation (3), 
Ri = — 


and 


P 2 i = (x c - D 2 )x + r { 
= ^Dx. + r 


( , i \ Pi , /*2 

( 7 + 1)^+7^ 


= -A[(7 J - !)/>* + ri 


Pi 


.a 


Pi + _P2_ 
Pli 3 P2i 3 


r, - 


x + r, 

- ^ r i] 

P2x 

(7 + 1 ) — 5 +7 — 5 
Pli 3 P2i 3 


Dx. 


In the rotating frame, r, and its time derivatives are given by 



Xj 


4'i - hi 


i I Vi ^ J V i J Xj 

r; = 

1 

1 

? r» = 

Vi + M 
2» 

, ri = 

Vi H~ 2J Xj J Vi. 

Zi 


Note that vectors p u , p 2l , and r, are dependent upon the eccentricity 
of Earth’s orbit, as is 1). To avoid this dependence, redefine these vectors 
relative to reference positions of Sun, Earth, and L 2 along the line of syzygy. 
Denoting the reference distances by an overbar. 

R, = x e x + r, ; 
u> c = fi 

Ri = he y + r* 

Ri = -f 2 X e X + he y + ft- 


Again using the definitions of 7 and p, x e — (7 + 1 — p)D. 
Using this form of the acceleration, 


— (7 + 1 — p)(~J x + Jy)D + r, 


Pi p 2 

Pli 3 P2i 3 


Vi - 


/ . -1 \ Pi . P2 

(7 + 1) — f + 7 — 3 , 
Pli 3 P2i J 


Dx, 


where 


Pii = (7 + l)Dx + r< and p 2i — jDx + r,. 


( 4 ) 
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Vectors r« and r, ; take the same forms as before, now relative to the reference 
L -2 location. 

For the relative motion, let 

r = r t - r,, = R, - R h . 


The 


a. 


r = -hi 


h 

Pir 




hi(7 + 1) ( 

Pit 


, T 

} 

J>2t 

1 


r h 


P'2/i 


[>\h 


1 1 

+ / ( >7 ( — 7 7 

P2/' 5 P2/C 


(5) 


Dx. 


as in Equation (1). 


2.3 Lunar Gravitational Effects 

This section discusses the lunar contribution to the telescope motion relative 
to the hub. Terms corresponding to the gravitational force of the moon upon 
the hub and telescope spacecraft are included in the relative equations ol 
motion (telescope relative to hub). The resulting contribution is then cast 
as an additive perturbation to the elliptical restricted three-body problem 
equations, such that when the lunar motion is ignored, the contribution re- 
verts to the lunar mass placed at the Earth position. The moon is treated 
as a point mass. 

As shown in Figure 3, let p :ii denote the vector from the moon to space- 
craft q either the hub or a telescope, and let /i 3 denote the lunar Kepler ian 
constant. The lunar force per unit mass on spacecraft i is then given by 


“/'3 


P:a 

P3i 3 ’ 


( 6 ) 


Therefore, the contribution to the equations of motion of the telescope rela- 


tive to the hub is 





Psk \ 
f >M ?) ’ 


■(7) 


where the subscripts h and t refer respectively to the hub and telescope. 
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2.4 Solar Radiation Pressure Effects 

The work presented in this section and in Appendix C was prepared by Pro- 
fessor David Richardson of the University of Cincinnati, under subcontract 
for this project [3]. 

The pressure from solar radiation imparts a tiny force on a telescope 
spacecraft. Depending upon spacecraft design and distance from the sun, 
the force can perturb both the spacecraft’s attitude and orbit. The force can 
also be harnessed to beneficially propel the spacecraft. 

We present a model to compute the force on a spacecraft, accounting for 
the force reduction when the spacecraft orbits through the terrestrial shadow. 

At a distance p\ t from the sun, the solar flux / (the irradiance) acting on 
the spacecraft is given by 

/-A. 

4np 2 u 

where 

Ij = 3.842 x 10 26 watts 
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is the solar luminosity (total omitted radiation). 

If A is the cross-sectional area of a spacecraft of mass m projected normal 
to the spacecraft-Smi line, then the solar radiation force F s per unit mass 
acting on the spacecraft is 


c r la 

Amncp 2 , 


1.0198 x 10 i7 CrA 
rn P 2 )i. 


N/unit mass. 


where c is the speed of light and 0 < Or < 2 is the parameter characteristic 
of the reflectivity of the spacecraft surface facing the sun: 


Cr = 0 translucent, 

Cr — 1 perfectly absorbent, 
Cr — 2 perfectly reflective. 


For trajectory motion that passes through any portion of Earth’s shadow, 
the full disk of the sun will be partially obscured. In the vicinity of L< 2 , 
this will occur at distances normal to the line of syzygy of approximately 
13,420 km or less. In such cases, the force expression above must be scaled 
by a “luminosity reduction factor 5 a which ranges from zero (total eclipse) 
to unity (full sunlight). The appropriate expression for the force per unit 
mass is then written: 


F. s = 


1.0198 x 10 17 C R Aa 
m Pu 


Pit = 


1.0198 x 10 l7 C R Aa 
w 1 1 (7 ~ l)Dx + iy|| 3 


(7+ l)Dx + r t _. 


( 8 ) 


The calculation of the luminosity reduction factor, a : is presented in Ap- 
pendix C. 

Our solar model did not consider the 11-year solar cycle, or estimate daily 
variations of solar flux, or the geomagnetic tail. 


2.5 Spacecraft Thruster Effects 

This section presents the derivation used to incorporate the effects of body- 
mounted thrusters in our equations of motion. 

First, we give the body-fixed acceleration components imparted by the 
thrusters. These terms are calculated as shown: 

x b = F x /m, jjb = Fy/rn, z b = F z /m. 
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A 

x b 
z b 

Figure 4: Vehicle Thrust in Body-Fixed Frame 

where F x , F y , F z are the components of the thrust F thrust in the body-fixed 
frame, as shown in Figure 4; m is the vehicle mass. 

Consider an arbitrary alignment of a body-fixed coordinate frame with 
respect to the rotating x, y , £ frame. The body-fixed yt,, Zb coordinate 
frame has its origin at the spacecraft’s mass center. 

The components of thrust expressed in the body frame must be trans- 
formed into the rotating reference frame in order to correctly incorporate 
these forces into the description of the motion. One way to express this 
transformation is as follows: 

X Xb 

y =T yb , 

_ ^ J [ _ 

where T is the transformation matrix formed as a combination of Euler ro- 
tations. 

It is somewhat easier to visualize the individual rotations by considering 
the inverse rotation description, rotating instead from the rotating reference 
frame to that fixed in the spacecraft: 
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Figure 5(a) shows the desired orientation of the body frame with respect 
to the rotating reference frame. However, by first assuming the coincident 
alignment of both the body frame and the rotating reference frame, we de- 
velop the inverse transformation matrix T ~ 1 as a combination of Euler rota- 
tions through a set of Euler angles: 

1. first rotate angle (V» + /) about the z, b axis as shown in Figure 5(b), 
where / = f(t - to) 

2. second rotate angle 0 about the new orientation of the iji, axis, as shown 
in Figure 5(c) 

3. third rotate angle <p about new orientation of the axis, as shown in 
Figure 5(d) 

Combine the sequence of rotations as follows with a right- to-left ordering 
of the rotation matrices: 



I 

o 

o 

t— 1 


~C(0) 0 -3(B)' 


c(i> + f) s (•(/’ + /) o 

T -1 = 

0 c(<p) s(<p) 


0 1 0 


-s(i’ + /) c:(v* + /) 0 


0 -sit) c(4)_ 


s(0) 0 c-(6>) 


0 0 J 


= T-\<p)T- 1 (6)T-\p + f), 


where the functions c and s represent cosine and sine, respectively; the ro- 
tation matrix T~ x refers to the required Euler rotation matrix about the 
4 -axis. Finally, the desired transformation matrix T (from the body-fixed 
frame to the rotating frame) is expressed as the inverse of T~ l : 

T=(T- 1 )- 1 = T. ( '</’ + f)T v {9)TM). 

These expressions premultiply the thrust force per unit mass acting upon the 
telescope spacecraft. 

Consider the following example that demonstrates comparable thrust and 
solar forces: A 500 kg mass telescope spacecraft with a 1 mN thruster can 
produce an acceleration of 15 km/day 2 . By choosing a solar flux reflectivity 
parameter Cm of 1.5 and placing various surface areas normal to the sun, the 
spaceemlVs acceleration due to solar radiation pressure is shown in Table 1 1. 

Contrast these low magnitudes to accelerations due to terrestrial and so- 
lar gravity (as given by Equation 5) of about 4,800 km/day 2 . Of course, these 
magnitudes depend upon example. One point is that both the low thrust of 
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(a) frame orien- (b) ib -f / about z\> 

tat ions 



(c) 0 about yt, 


(d) (p about Xfr 


Figure 5: Frame Rotation from Body to Rotating Frames 


electric propulsion and the low pressure of solar radiation produce acceler- 
ations very much lower than the gravitational effects. Another important 
point is that solar radiation pressure can be used for orbit control because 
the force from solar radiation pressure may be comparable to that from a 
thruster suite. 

2.6 Summary 

In this section, the force equations derived in Sections 2.1 through 2.5 are 
combined. The resulting equation models the elliptical restricted three-body 
problem, incorporating lunar, solar radiation, and thrust perturbations. 
Combining the expressions of this section, given by Equations (4), (6), 


121 


Table 1 : Acceleration due to Solar Radiation Pressure vs. Surface Area 


Sun-Facing 

Solar Radiation Pressure 

Surface Area 

Acceleration 

(nr) 

(km/ day 2 ) 

1 (« 3 fi x 3 ft) 

1 

100 (« 33 ft x 33 ft) 

10 

150 (« 40 ft x 40 ft) 

15 

1,000 (wl 00 ft x 100 ft) 

100 


and ( 8 ), along with the thrust expressions of Section 2.5, the differential 
equations of motion of the telescope are given by: 


Rt = (7 + 1 - p){~f 2 x + fy)D - r ( 


(M Pi 

Pit* P2t* 

Pit 


r t 


Pit 


Pit* 1 


Dx 


Pi j 
Pit 

1.0198 x 10 l7 C R Aa 


m 


+ l)Dx + r,| 


3 [(7 + l)Dx + r t ] 



W) 

-m 

0 “ 


' c(0) 0 

«{Q) 


T o 

0 

+ 

m 

cm 

0 


0 1 

0 


o c(d>) 

~ s i4>) 


0 

0 

1 1 


-M0) 0 

r(0) 


0 .s(0) 

((&)_ 


thrust 

m 

where yj f = P + /. The corresponding equations for the hub are 
R/,. = (7 + 1 - p)(-/ 2 x + fy)D + r h 

a, a.-, 1 _ 

Dx 


thrust 



y , i \ p i , p-2 ' 

( / + 1 ) 3 + / o 

\P\h- Pih J 

Plh Plh _ 


P Ih 

Pi 7 

Plli 

1.0198 x IQ 17 C R Aa 
mj|( 7 + l)Dx + r/,| 


3 [(7 + l)-Dx + r fc ]. 


ER.3B 

lunar 

SRP 


Combining the relative motion expressions of this section (Equations (5) 
and (7)), along with the solar radiation pressure of Equation ( 8 ) and the 
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thrust expressions of Section 2.5, the differential equations of relative motion 
for a telescope spacecraft are given by: 


( r 4 v h 

r = -Ati — 3 5 

VPlt Plh 


— /X'2 


r h 


Pit Pih 


- Pi 


, w 1 1 

A*i (7 + 1) — 3 ; 

Pu s P\h- 


Pit Pih 


+ A'27 — o - 

Pit 3 


Pih 


+ 


,P3£ Pih. 
1.0198 X 10 17 CrA<j 


m||(7 + l)Dx + r t 
1.0198 x 10 17 C R Aa 


3 [(7 + l)Dx + r f ] 


m||(7 + l)Dx + rfc|| 


3 [(7 + 1)-Dx + r/,] 


Dx 



'cU>') 

-s (i’ 1 ) 

0" 


c(d) 

0 

■m] 


"l 

0 

0 

+ 

» W) 

r(V') 

0 


0 

1 

0 


0 

r-W 

-•*#) 


0 

0 

1 


-*(8) 

0 

c(0) 


0 

s(4>) 

m. 


thrust 
777 , 


ER3B 

lunar 

SRP f 

SRP fc 

thrust 


We present the following example to permit comparison of the relative 
contribution of the terms to the telescope motion. The initial conditions are 
listed in Table 2: they were taken from the examples discussed in Part 2 of 
the research. As mentioned in that report, these values were selected so as 
to excite only the oscillatory linear modes. 

Figure 6 presents the solution to numerical integration of four different 
force models selected from the summary equation above. The models were 
applied to both the hub and telescope, and the resulting state vectors were 
differenced in order to determine the position of the telescope relative to the 
hub. In each case, the same force model was applied to both the hub and 
telescope. Additionally, the same value of the reference distance D was used 
for all force models. 

The reference solution is represented in the figure by the x-axis. This 
refers to the circular restricted model case. The circular restricted model so- 
lution is obtained by including the elliptical restricted model (ER3B) terms, 
but with Earth’s orbital eccentricity treated as being zero; this duplicates 
the model discussed in Part 2. The other models depicted in Figure 6 rep- 
resent the addition of ellipticity (ER3B), lunar, anti solar radiation (SRP) 
effects over 20 days. It is noted that the elliptical contribution provides the 
dominant perturbation to the circular restricted solution. In this example, 
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Table 2: Example Initial Conditions 



hub 

(state rel. 
to L>) 

telescope 
(state rel. 
to L>) 

telescope 
(state rel. 
to hub) 

,r(0) (km) 

-227,219.419 

-227,219.483 

-0.064780 

£/(0) (km) 

0.0 

0.0 

0.0 

2 ( 0 ) (kin) 

-250,000.000 

-249,999.974 

0.026445 

.r(0) (km/day) 

0.0 

0.0 

0.0 

i/(0) (km/day) 

25,625.039 

25,625.044 

0.004421 

£(0) (km/day) 

0.0 

0.0 

0.0 

mass m(kg) 

500 

500 


F tln-iist (N) 

0 

0 


sun-facing 
area A (in 2 ) 

150 

150 



should the viewing of a scientific target be limited to no more than 10 days, 
the perturbation due to the elliptical contribution is less than 1 nr. 

In the MAXIM or similar missions, there may or may not be an actual 
hub spacecraft located at the aperture's center. Regardless, it is necessary 
to treat, the hub as a central reference point for locating the positions of 
the individual telescope spacecraft. However, unlike the gravitational forces 
present, the solar radiation pressure effect upon a spacecraft with actual 
mass and area is substantial as compared to that upon a. ‘phantom" hub. 
Therefore, for purposes of the simulation, the hub “spacecraft is treated as 
having the same physical characteristics as the telescope. In this manner, the 
hub is maintained as an adequate reference for the position of the telescope. 

The effects of thrusters were not simulated here, due to the vast uncer- 
tainties of force (both magnitude and direction) and duration. We considered 
trade studies involving thruster forces to be beyond the focus of this investi- 
gation. 
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Circ Circ+Ellip Circ+Ellip+Lunar Circ+Ellip+Lunar+Solar 


Figure 6: Effects of Perturbations on Relative Distance — Full Equations 
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3 Expanded Equations of Motion 

3.1 Circular Restricted Three-Body Problem 

In Part 2 [2], Equation (1) is expanded through terms which are linear in 
the coordinates of r and no more than cubic in the coordinates of r/,.. This 
expansion takes the form 


r = A :i [— r + 3xx] 

+ A 4 [3xr h . + 3x h r + (3r h ■ r - I5 .x.t,,.)x] 

+ A s [(3r fc • r - 15.vx h )r h + |(r h 2 - 5x h 2 )r 
- y(2.x, l r, 1 ■ r - 7xx h 2 + xr h 2 )x ] . 

where the constants . 1 7 are given by 

a.. = ^ ^ 

(x e + Di)' (x e — D 2 Y 

Terms involving A z are considered to be of order z; terms of order lower than 
3 do not appear. 

The acceleration vector r may be written relative to a rotating coordinate 
system which rotates at the constant angular rate n about the z-axis normal 
to the ecliptic, and with the x direction as previously defined. This gives 


r — 


x — 2m/ — n 2 x 
y + 2 mi — n 2 y 
z 


where the column vector notation is used to indicate the xyz vector compo- 
nents. 


3.2 Elliptical Restricted Three-Body Problem 

Consider the elliptical restricted three-body equations of motion as given in 
Equation (5). The right side of this set of relative acceleration equations may 
be expanded as in Part. 2 [2]. Consider the effects of various contributions to 


126 



the magnitude ordering scheme. For ordering purposes, take 


(r h == 600. 000 km) 

(r = 0.5 km). 

A rough estimate may be obtained of the contributions that the various 
perturbations make to the circular restricted problem solution. Say that a 
perturbing term may be treated as modifying the linear frequencies asso- 
ciated with the circular restricted problem, and consider the square of the 
perturbing frequency to be roughly the magnitude of the coefficient of r in 
the perturbing acceleration. (Recall from Part 1 that the linear periods in 
and out of the xy-plane are approximately 177.566 days and 184.002 clays, 
respectively.) Then, after 90 days, the effects of the terms containing various 
powers of r and r h are given in Table 3. with effects included of roughly 20 m 
and larger. 

Table 3: Along- Ellipsoid Effect of Sun-Earth Perturbation on Solution (90 
days) 


perturbing 

term 

effect 

(m) 

re 

16.7 

rr h 

289.3 

rr h 2 

118.5 

rr h 3 

47.4 

rr h 4 

18.9 


p = 3.04 x 10~ 6 
7 = 1.01 x 10“ 2 

^ = 4.0 x 10" a 
= 3.2 x 1CT 9 


Wo recommend retaining terms which contribute down to approximately 
20 rn (can keep additional /q terms if convenient). Retaining the rr^ s terms 
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results in the expanded/tmnrated differential equations 


-f A 4 [3-rr/j + 3x h r + (3r h • r - 15xx h )±] 

+ A- a [(3r h • r - 15xx h )r h - fO’/. 2 - 5x h 2 )r 
- ] -§(2x k r k ■ r - 7xx h 2 + xr h 2 )x] 

+ A 6 [fr h{-xr h 2 - 2x h r ■ r h + 7xx h 2 ) + |r(7ar & 3 - 3 x h r h 2 ) 
+ y(-r h 2 r A ■ r + 7xr h 2 x h + 7x h 2 r ■ r h - 21;r.rp)x] . 


Next consider the left side of the differential equation: 


r = 


;• - fy - 2 iii - f 2 x 

jj + fx + 2 J'x - J' 2 (j 


— 2 n c y — n c 2 x 


-fy -2 (/ - n c )y 

- (/ 2 - n 2 )x - 

+ 2n r x — n, 2 y 

+ 

fx + 2 (/ - n c )x 

nc 2 )y 

z 


0 



( 10 ) 


In tins expanded form, the first vector term represents the acceleration which 
appears in the circular restricted problem (n c refers to the mean motion of 
the circular restricted Earth orbit). The second vector term gives the per- 
turbation which is added by including the elliptic restricted effects. The 
perturbation term is to be expanded in terms of the eccentricity e of Earth's 
orbit (« 0.017). In keeping with the earlier magnitude ordering, it is esti- 
mated that it is sufficient to retain only contributions which are linear in e, 
because er is approximately 8.5 m. 

Say that the circular restricted problem takes D as bei ng the reference 
value D. Then, the corresponding mean motion is n (: = vp//3 a . Now, in 
the elliptic restricted problem, the mean motion is 


n — 




Ke- 


if IJ is chosen to be the semi-major axis u, then n = n (: . If D ~ D m (;an = 
ci( 1 -f- e 2 /2), then 
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Expressing / and / in terms of time. 

• h v///,a(l — e 2 ) na 2 \/l — e 2 /a\ 2 

/=D5 = ^ =— D —« n A-d) ■ 

where n — as given by Kepler : s third law. Using the two-bodv 

relationship 

D = q(l - e 2 ) 

1 + fi cos / 

• / 1 + ecos/\ 2 . 

/ ~ «c ( 1 _ e2 J ~ n c (l + 2ecos /). 

In terms of the mean anomaly £, 

f « n c (l -I- 2ecos£). 


Differentiating, 


/ as -2e?r c 2 sin t. 


Using these results in the perturbing vector of Equation (10) gives 


-fy - 2(J - n c )y - (/ 2 - n*)x 
!■<■ + 2 (./' - n c )x - (f 2 - n 2 )y 
0 



—2< j yn,. 2 sin — Aeyn,. cos C. — \< j .i:n, 2 cos/' 
2exn c 2 sin £ — 4ex?i c cos (‘ — 4 eyn c 2 cos / 


0 


3.3 Lunar Gravitational Effects 

Consider the lunar force given above by Equation (7). 

As in the previous work, it is desired that this contribution be expressed 
in terms of the position of the telescope relative to the hub. Let the vector 
r again refer to the relative telescope position. Therefore, 

Put. = Pzk + r - ( n ) 

As in the earlier analysis of the effects of Sun and Earth, the contribution of 
Equation (7) may be written as an expansion in terms of r and its compo- 
nents. Accordingly, a similar binomial expansion development is followed as 
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that employed in Part 2 of the research. The square of the magnitude of p u 
is given by 


Then, 


l>:n 2 — P:n ' Pit 

= Pih 2 + r ' 2 + 2 Pih ■ r 



2 Pih • r 

Pih 2 


Pit" 


Pih' 


1 + 


Pih 


-3 


(- 

\pm 

(i + «" V2 


+ 


-Pih ■ r 

Pih 2 


1 - 3/2 


where 




2Pa/, ■ r 

Pil 2 


is assumed to be less than unity. Using a. binomial expansion, 





Using this expansion in Equation (7), and using Equation (11) to substitute 
for p M gives 


~ Pi' 


Pih" 


> + E 


fc=l 


-3/2 V, 
k ^ 



Expanding through linear terms in r, 


F 

m 



+ 3 /<3 


Pih r ~ Pih 

pill 1 Ps/i 2 


( 12 ) 


It is preferable to treat these terms as an additive perturbation to the 
equations of the elliptical restricted three-body problem. Assume that the 
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baseline elliptical restricted problem contains an object with combined terres- 
trial and lunar mass, located at the mean Earth-Moon barycenter. Therefore, 
the expansion of Equation (7) should contain the lunar contribution to this 
combined mass along with terms representing the effect of the lunar motion 
about the barycenter. Accordingly, Equation (12) may be rewritten as 


~Vi 



3 Pih r ' Pih \ 
Pih 2 P'lh 2 J 


~ Pi r 



1 

P2li’ J 


+ 3p 3 r • 


( PihPih 

Pik 5 




PlhPlh\ 

Pih 5 J 


(13) 


Because the vector p 2 /< was used in the earlier analysis, it is convenient 
to write p ih as 

Pik = Pih + r,. m , 

where r em refers to the lunar position relative to the mean barycenter. The 
square of the magnitude of p ih is given by 


Pih Pih * Pih 

— P’2h~ + J'em 2 + ^Plh 



+ 


^Pi)) • r £ 

Pih 2 


This gives 


J_ _ J_ 
Pih P'lh 


1 + 


J cm 

P2h 


9 


*+ 


P2h ’ 

Pih' 2 


- 1 /2 


Using a Legendre polynomial expansion. 


— = -Lff- ] n-( coss), 

P3k P2h “J V P2A 


where 


c A P‘2h * r e?n 

cos 6 = . 

P2lp'ern 

Expanding through terms linear in r cm . 


_L ~ 1 

P'Sk P2h 


— ( 1 + —■ cos S ) , 

Pih 
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giving 


P'Sh 

1 


3 


p2h 


1 


, o' em C r 

1 + 3 — COS O 

p2h 

i , r ? em C* 

1 + 5 cos 5 


P5h* P2h° \ P2h 
Using these representations in Equation (13), 

, r 3p 2h r • p 2/l 

~Pi 


P'2li P2h' J 

3 P'3 [^P'2/i ' *"em “I” 1 ’ P 2 h^Gm Y ’ 1 emp 2 h 
— 5l’ • p2hP‘2hP2h ‘ r e?n/ P2/i ] / p2h~ • 


(14) 


Examining Equation (14), the form of the first term in F m is identical 
to the corresponding terrestrial term derived in the earlier work. The only 
difference is that, here, the mass coefficient is /z ;i rather than /z This term 
represents the effect of a lunar mass collocated with the terrestrial mass; 
the remainder of F m represents the effect of the lunar motion about Earth. 
The first term may be added to the earlier work simply by replacing //. 2 
with (i 2 + /y in the relative equations of motion for the elliptical restricted 
three-body problem. 

For the second term of Equation (14), the Earth-hub position vector p 2h 
is expanded in a fashion similar to that of the earlier work. In the context of 
the elliptical restricted problem using a reference location of L 2 , 


P2h = 7 Dx + r h 


where 7 and D are constants as previously defined. The square of the mag- 
nitude of p 2h is given by 


P2h — P-2h ' P‘2h 

, __ 2 


— (7 D)~ + + 2"/ Dr.h 


= "fir 


1 + 'P5 


2 Xh 

ID 


where xj h is the x-component of r /; .. This gives 


1 /p 2 h = (1 + r h ) 1 ^/( 1 D). 


where 


2 


A 


£h = 



2 Xh 

I'D 


Once again, powers of this fraction are formed using binomial expansions, 
giving 



The negative powers of pok that appear in the second term of Equa- 
tion (14) are then formed using the appropriate values of l in this binomial 
expansion. Expanding through linear terms in r h and substituting, F m be- 
comes 


“/* 3 


3p2/i r ' P'2h 


„P2hr p2h 

3 P'S [( 2x.T em 4“ yy em 4^. e7H )x -(- ( yX-cjn T y 

+ + X2 em )z] /(7D) 4 

- yu 3 [- 15x em x /l r + 3r,, • r em r - 15xr em x h 


+ 3r ■ r fe r,. Ttl + 3(r • r ( , m )(-5x,,x + r h ) + lObxx^XhX 
- lFjxx ( . m r h - 15xr* ■ r cm x - 15,r m r • r k x]/(jDf. 


(15) 


Again, a rough estimate of the contributions that the various perturba- 
tions make to the solution are presented. As for the inclusion of the ellipticity 
in Section 3.2, the perturbing terms are treated as modifying the linear fre- 
quencies. After 90 flays, the effects of the lunar terms containing various 
powers of the relevant variables are given in Table 4, including contributions 
greater than roughly 0.9 111 . 


Table 4: Along-Ellipsoicl Effect of Lunar Perturbation on Solution (90 days) 


perturbing 

effect 

term 

(m) 

? ^ em 

2.3 

rr 2 

0.6 

1 1 h? em 

0.9 
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3.4 Summary 

In this section, the force equations derived in Sections 3.1 through 3.3 are 
combined. The models for solar radiation pressure and thrust are repeated 
here from the equation of Section 2.6 for completeness. The resulting equa- 
tion models the expanded elliptical restricted three-body problem, incorpo- 
rating lunar, solar radiation, and thrust perturbations. 

Combining the expressions of this section, the expanded differential equa- 
tions of relative motion for a telescope spacecraft a, re given by: 


r = /l 3 [-r + 3xx] + A, { [3xr/, + 3x h r + (3r/ t • r - 15x.r*)x] 

+ .4 5 [(3r h ■ r - 15xx h )r h + \{r h 2 - 5x h 2 )r 
- f (2x h v h • r - 7xx h 2 + xr,, 2 )x] 

+ /In .w/, 2 - 2x/,r • r h 4- 7x.v h 2 ) + §r(7x fc 3 - 3 x h r h 2 ) 
+ y (~r h 2 r h ■ r + 7 xr h 2 x h + 7x, ( 2 r • r* - 21xxp)x] 

3/^3 [( 2 xx ern + y y cm 4 " »*! bii)x 4 ” ( yx em 4 " xy ein )y 

+ {zx ern + 2* f)n )zj / (7 D) 

[ 15x CTrt .r/ ( r 4" or/ t ■ r CTn Y loxr cm .x/j 

+ 3r • r,,r rm 4- 3(r • r, !m )(-5x,,x 4- 17 ,) 4- 105xx,.„,x,,x 
- 15xx em r,, - 15.ix h • r em x - 15x em r • r h x]/(~/ Df 


+ 


x.uiyo a iu 


7 + l)Dx 4- r,, 


m||(7 4- L)L>x + r>, 4- r| 

1.0198 x IQ 17 CrA(t r . , 

HI(-, + l)5x + r,.||4 (-> + 1)Cx + r -] 
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where y/ = ip + /. 

We present the following example to permit comparison of the relative 
contribution of the terms to the telescope motion. The initial conditions are 
the same as those presented earlier in Table' 2. 

As in Figure 6, Figure 7 presents the solution to four different force mod- 
els selected from the relative motion summary equation above. Where hub 
position was required, it was obtained from separate integration of the full, 
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Circ Circ+Ellip Circ+Ellip+Lunar ~ - Circ+Ellip+Lunar+Solar 


Figure 7: Effects of Perturbations on Relative Distance — Expanded Equa- 
tions 

unexpanded hub equations from Section 2.6, using the same force model as 
for the relative motion. 

Once again, the reference solution is represented in the figure by the x- 
axis. This refers to the circular restricted model case. The other models 
depicted in Figure 6 represent the addition of ellipticity (ER3B), lunar, and 
solar radiation (SRP) effects over 20 days. Solar radiation is again treated as 
being applied to both the telescope and to either a real or phantom spacecraft 
at the hub, with the same physical characteristics as the telescope. 
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4 Modeling Uncertainties 

4.1 Elliptical Restricted Three-Body Problem 

Very Similar Values for Different Definitions. In computing the value 
of /), the reference distance between the centers of the sun and earth, we need 
to determine the mean distance between Earth's center and the Earth— Moon 
bary center. 

The astronomical unit (AU) is defined in Scidelmann [4] as the radius of 
a circular orbit in which a body of negligible mass, and free of perturbations, 
would revolve around the sun in 2 ir/k days, where k is the Gaussian grav- 
itational constant. This is slightly less than the semi-major axis of Earth’s 
orbit. Begin with the distance of the AU: 

AU = 149,597,870.000 km [4, page 700, IAU System] 

— 149,597,870.660 km [4, page 700, Best Estimate] 

Yet compare these values with that from Dunham and Muhonen [5, 
page 200]: 

AU - 149,597,870.691 km 

The mean distance from Earth to the sun also has different values: 

1.0000010178 AU [4, page 700, IAU System] 

1.00000105726665 AU [4. page 700, Best Estimate] 

The calculation of the actual mean distance yields the following: 

149, 598, 022.261 km, using the values from the IAU System 
149,598,028.825 km, using the values from the Best Estimate 

Again, compare these values with that from Dunham and Muhonen [5, 
page 200]: 

Earth+Moon semi-major axis = 1.000001018 AU 

= 149. 598. 023.0' km 

Dunham and Muhonen define their value as the mean distance of the 
Earth+Moon barycenter from the sun. Contrast this value to that of Seidel- 
mann, who specifies this value as the mean distance of only Earth from the 
sun; however, the values are very close to one another. 


Wo continue using those prominent resources for the eccentricity, e. 

On page 700, Seidelmann gives the value of 0.016708617 for the mean 
eccentricity of Earth's orbit about the sun. On page 200, Dunham and 
Muhoncn give the value of 0.01670862 as the Earth+Moon’s eccentricity 
about the sun. 

Sensitivity of D. From the above, the derived constant D depends on 
two constants with different values: 

D = a( 1 + e 2 /2) 

= 149,618,904.49 km, using IAU System for a and Seidelmann for e 
= 149,618,911.06 km, using Best Estimate for a and Seidelmann for e 
= 149,618,905.22 km, using Dunham and Muhoncn for both a and e 

These different values were applied in separate simulations and do not 
substantially affect the analysis results. 

4.2 Lunar Gravitational Effects 

Description of the Lunar Motion Around Earth. The moon’s motion 
around Earth is not specified by an analytical model. The motion is specified 
by the software and cphemericles files provided by the JPL. For the needs of 
this report, the JPL ephemerides provide the positions of Sun, Earth, and 
Moon to very high precision. The JPL ephemerides are given as blocks of 
Chebyshev coefficients, which, when interpolated, reproduce the original JPL 
numerical integrations to within 1.5 cm. 

The instructions for using the JPLEPH.200 ephemeris files state that one 
calls the JPL model with the moon as the target and Earth+Moon barycenter 
as reference point. Yet a comment was found within the FORTRAN code 
that the lunar state is always geocentric. 
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5 Sensitivity Summary 

It is recognized that the hub position knowledge is likely to be much less pre- 
cise than that of the relative motion, perhaps on the order of 1 km. There- 
fore, it is useful to consider the effect of this imprecision upon the telescope s 
position relative to the hub. 

Appendix A presents a detailed derivation of the variational equations 
used to address the sensitivity of telescope motion to errors in knowledge of 
the hub position. 

The results of the analysis clearly show that the errors in telescope posi- 
tion relative to the hub. based on knowledge of hub position to 1.7 km are 
sufficiently small that they may be ignored. 

Several sample tests of this behavior were conducted. For one test case, 
the initial conditions of Table 2 were used hero as a nominal set of initial 
conditions for the hub and telescope. 

The integration was then performed with the hub offset from its nominal 
initial state by 1 km in each directional component. This was to simulate 
an initial error in hub position. As seen in Figure 8. the telescope position 
relative to the hub differs from the nominal case by millimeters over 40 days. 

The simulation was also repeated with the initial hub position offset from 
the nominal state by 17.3 km (10 km in each direction). Figure 9 demon- 
strates the roughly ten-fold increase in error over the previous example. 
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Figure 8: Relative Motion Error Caused by 1.7 km Initial Hub Position Error 

6 Summary and Conclusion 

This report details our further work describing the formation flying between 
spacecraft near the Sun-Earth L 2 libration point, beginning with the circular 
restricted three-bocly problem for the hub motion about L 2 . 

Continuing from our previous works, these analyses develop the elliptical 
restricted three-body problem from previous work with circular problem. We 
built on our familiarity with the circular problem to address the following as 
perturbations upon the circular problem: 

• elliptical orbit of Earth-Moon about Sun 

• lunar gravitational effects 

• solar radiation pressure effects 

• thrusters on vehicle 

These were incorporated as additive perturbations to the circular re- 
stricted three-body problem with expansions of varying levels of fidelity. We 
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Figure 9: Relative Motion Error Caused by 17.3 km Initial Hub Position 
Error 


develop two forms of models of the full derivation: 

• full nonlinear baseline 

• expanded form of full nonlinear baseline to appropriate order 

Section 2 presented the derivation of the full nonlinear baseline with the 
perturbations. The equations were implemented through their coding in a 
MATLAB simulation. One example was presented using the discovery from 
Part 2 of valid initial conditions that excite only oscillatory motion in the 
linear modes. The results shown in Figure 6 demonstrate the dominant 
perturbation due to the elliptical motion of Earth about the sun. In this 
example, the perturbations due to lunar and solar radiation pressure are 
negligible. Overall, the figure indicates that during the typical 5-dav science 
observation of this example, these perturbations contribute less than 1 m 
of relative position difference. The proposed electric propulsion should have 
little trouble maintaining position. 


140 


Section 3 presented the expanded circular and elliptical portions of the 
full model through terms which are linear in the coordinates of the telescope 
position relative to the hub and no more than cubic in the coordinates of 
the hub position relative to /v 2 . The derivation describes the magnitude of 
the various terms and why some were truncated. Additionally, the lunar 
gravity model was expanded and some terms truncated due to very small 
perturbations upon the motion. As before, the perturbations due to the 
elliptical motion of the Earth about the Sun are dominant. Again, lunar and 
solar radiation pressure effects are negligible. 

Modeling uncertainties were discussed in Section 4. One type of uncer- 
tainty is that the same listed number can be found to have different defini- 
tions between two popular references. Calculations, of course, give slightly 
different results; however, they do not substantially affect the results. 

Section 5 explains the results of a longer derivation given in Appendix A. 
During the 18-months of part-time work on this report, our sponsor requested 
that we investigate the sensitivity of telescope motion' to errors in hub posi- 
tion. The analysis and simulation show that errors in the knowledge of hub 
position of 1 km and 10 km (in each of the three position components) yield 
millimeters of error in telescope position. 

Appendix B presents a back-of-the-envelope look at the relative motion 
between a hub and telescope at the Earth-Moon L2 point. Our sponsor made 
a comment that future work may be redirected to this vicinity. Our work 
addressed the effects of the largest contributions to the elliptical restricted 
problem. The perturbation caused by solar gravitation was found to be 
negligible. 

Appendix C presents a derivation of the luminosity reduction factor in 
the context of the solar radiation pressure model. The luminosity reduction 
factor accounts for shadowing effects due to a partial eclipse. 


In conclusion, based on earlier literature searches, we believe this new 
work is unique because it describes the primary perturbations to the de- 
scription of relative motion between nearby spacecraft. We verified that the 
effects of the elliptical motion of the Earth about the Sun is the dominant 
perturbation to the circular restricted three body problem. Contributions 
due to lunar gravity and solar radiation pressure are nearly negligible in the 
chosen examples. 
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Appendices 

A Sensitivity of Telescope Motion to Errors 
in Hub Position 

c This appendix demonstrates the low sensitivity to errors in hub position, of 

the telescope motion relative to the hub. The effects of 1 km hub position 
errors upon the relative telescope motion are presented. The model used is 
the elliptical restricted model as presented in Section 3.2. 

A.l Derivation of Variational Equations 

Consider the second-order differential equation of motion of a telescope rel- 
ative to the hub, assumed to take the form of the earlier expanded and 
truncated form: 

r = A 3 [— r + 3;rx] 

-I- A 4 [3-rrft + 3x h r 4- (3r h • r - I5xx h )x] 

+ A 5 [(3r, x ■ r - I5xx h )r h + § (r h 2 - 5x h 2 )r 
- ^(2x h r k • r - 7 xx h 2 + xr h 2 )x] 

+ An r h (-xr h 2 - 2x h r ■ r h + 7xx h 2 ) + |r(7.r^ 3 - 3;r h r h 2 ) 

+ ^{-r h 2 r h • r + 7 xr h 2 x h + 7x h 2 r ■ r h - 21x.t /i 3 )x] , 

where 

i 

r = [x y z \ 



x — 2 n c y — n c 2 x 


2 eyrie 2 sin £ ~ 4 eyn c cos £ — 4 exn c 2 cos l 

r = 

V + 2 n c x - n c 2 y 

+ 

—2exn c 2 sin f: + 4 e. xn c cos £ — 4 eyrie 2 cos £ 


z 


0 


A = — + **1 

(1 + 7 ) n D n 7 n D n ' 
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and 


r = telescope position relative to the hub 
r/, = hub position relative to L-z 
x, y, z — Cartesian components of r 
//] = solar Keplerian constant 

(.i 2 = terrestrial Keplerian constant (Earth + Moon) 

D — reference Sun-Earth distance 
x ( , = distance from system barycenter to L 2 
D i = distance from system barycenter to Sun 
D 2 ~ distance from system barycenter to Earth-Moon barycenter 
7 = (x e - D 2 )/D 

n c ~ circular restricted terrestrial mean motion 
t — mean anomaly of Earth orbit. 

The rotating coordinate system is defined to have .r-axis along the Sun-Earth 
line of syzygy and z-axis normal to Earth's orbit about the sun. 

The vector differential equation may be written in the form 


:r 

y 

= f = 

■ /,: ‘ 
fy 



fz 


The vector function f represents the acceleration relative to the rotating 
frame. This function is linear in x, y, z, and their derivatives. 

Say that there is a nominal solution r„ associated with a nominal hub 
position ry n . It is desired to estimate the sensitivity of the relative telescope 
motion to errors in r/,. which are on the order of 1 km. 

Define the state vector x as 

x = [ x y z x ij z ] 1 . 

The time derivative of this state is then 

X = [ X y i f x fy /- ] 1 . 

Referring to this vector as the function F(x.r/,). the differential equation 
now takes the first-order form 

x = F. (16) 
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where the nominal solution satisfies 


= F 

— x n 


(17) 


The solution x is now considered to be x n as perturbed by an error in 
Th. It is desired to examine the effect of this error upon x. Writing x as 
a linearized Maclaurin series in this error (or, equivalently, a Taylor series 
about the nominal hub trajectory) gives 


where 
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( 18 ) 


The subscript n refers to evaluation using the nominal trajectories; Xh, yn, 
and Zh arc the Cartesian components of the hub motion relative to £ 2 - 

It is the quantities of the elements of this matrix U that are of interest, 
in particular, the partial derivatives of the position components. Once these 
elements are determined, the effect of the hub position error upon the relative 
telescope motion may be approximated by 
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Because the components of the hub position error are presumed to bo no 
greater than 1 kin, it is sufficient to examine the position partials within 
U ; they are treated as being scaled by unity to form the relative position 
variations. 

A Maclaurin series may be formed as well for the derivative of the state. 


as 


x = x n + 


dx 

Ovi, 


{j-h I*/m) * 


This derivative may instead be formed by differentiation of Equation (18), 
giving 




* ( r /i — r/m) + 


dx 

dr h 


dt 


( r - T hn ) 


(19) 


It is assumed that the change in hub position error is of higher order than 
the error itself. Therefore, these two equations 'indie -ate that 


dx 

dr h 


d { dx \ 
dt V dr h J ’ 


( 20 ) 


The function F may also be written as a linearized series, giving 


F as F„ 


OF 

. . dF 

dx 

0r h 

■ [r h - Thn) + 7T- 
n (JX 

n 0r >’ 


(l /?. 1 hn) : 


( 21 ) 


recognizing the dependence of F upon r* both directly as well as indirectly 
through x. Substituting Equations (19) and (21) into Equation (16), and 
taking into account the nominal solution of Equation (17) as well as the 
approximation of Equation (20), 


d I 

( c)x 

\ OF 

dx 

OF 

1 

dt ' 

\ Or i, 

J “ Ox 

» 0r h \ 

dr/, 

n 11 


or 




( / + 2 £ 

or /( 


( 22 ) 


Recall that the vector function F is linear in x, with no additive constant. 
Therefore, 

F = — • x 

dx 
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Substituting into Equation (17) 


Xn 


dF 


dx 


■ x n . 

n 


Note the similarity of this equation to the homogeneous part of Equation (22). 
This indicates that any column of the homogeneous solution of Equation (22) 
takes the form of the general solution of Equation (17). 


A. 2 Linear Equations 


It may be seen that, for the purpose of determining the gross effects of hub 
position error on the relative telescope motion, it is sufficient to consider 
only the first-order contributions to the differential equations. Say that the 
nominal solution for x is known in the form of a Taylor series. Therefore, 
the homogeneous solution to Equation (22) is also known in the form of a 
Taylor series, given as 

e 2 

U h — Uk o + £ Uhi + 7T^2 + • • • • 

Similarly, let the nominal system matrix also be written as a Taylor series: 


d F 

! hi 


dF 

dx 


+ e 


nO 


dF 

dx 


n 1 


SiTF 

+ Hdx 


+ 


77.2 


Because the dominant terms of F are void of the components of rj,, the 
nonhomogeneous part of Equation (22) is of higher order, with Taylor series 


dF 

dF 

e 2 dF 

dr h 

n ~ f < )r h 

nl + "2 foh 


Now, consider the particular solution to Equation (22) as taking the form of 
a Taylor series as well. Because the nonhomogeneous part of the equation 
begins at order 1, so does the particular solution: 

U p ~ eUpi + — U P 2 + • • • . 


Substituting these series expansions into Equation (22) gives, at order 



As previously mentioned, this solution is assumed to be known from the 
knowledge of the nominal motion. At order one, 



• U h i + 

tlO 


OF 


<9x 


• U h0: 
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the solution to which is again assumed to bo known, and 


Ur ' dx 


n 0 


„ OF 
• C'l d" 77 


n 1 


(24) 


For Equation (24), only the particular solution is required. The system ma- 
trix here is the same as that of Equation (23); terms of the homogeneous solu- 
tion to Equation (24) are already provided in the solution to Equation (23). 
It is assumed that the solutions are convergent for any given time; there- 
fore, the gross contribution to the particular solution may be found from this 
equation alone. Accordingly, henceforth, the matrix U is assumed to refer to 
those terms through LA, j and f/ 7 u. 

For purposes of this analysis, the order 0 terms are considered to be those 
with the numerical coefficient A 3 ; order 1 terms are those with the coefficient 
A] . Other terms, including those involving the eccentricity e are treated as 
higher order. For that reason, the mean motion n of Earths orbit may be 
approximated by n c . 


A. 3 Solution Development 

The solution for each vector of the linear system is constructed in the stan- 
dard fashion. Let u represent a column vector of U , and v represent, the cor- 
responding column of the nonhomogeneous part of Equation (22). through 
the prescribed order. Then, the linear system takes the form 

u = Au + v. 

where the matrix A refers to the matrix of partial derivatives seen in Equa- 
tion (23). The homogeneous solution may be constructed using standard 
methods, as will be seen below; therefore, it is necessary to develop the par- 
ticular solution to the system of differential equations given by Equation (24). 

In these differential equations, the nonhomogeneous part contains trigono- 
metric terms, which may be written as exponential terms with exponential 
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multipliers equal to the characteristic multipliers of the homogeneous system. 
Accordingly, the particular solution is developed for the case of a single such 
exponential term at a time; the solution for a number of such terms may be 
constructed by superposition of the individual solutions. It is assumed that 
the characteristic multipliers are distinct, as is the case for the hub-telescope 
system. 

A. 3.1 Non-secular Solution 

First, examine the case where the nonhomogeneous vector v contains periodic 
(or, equivalently, exponential) contributions with frequencies not equal to one 
of the eigenvectors of the system matrix A . 

Consider a general first-order nonhomogeneous linear system of differen- 
tial equations given by 

x = dx + e lSlt ej, 

where e 0 is the j th system unit vector, and id is not one of the system 
eigenvalues. The particular solution of such a system may be found using 
the method of undetermined coefficients. Say that the particular solution is 
of the form 

x p = be tClt . 

where the vector b is a constant vector which is to be determined. Substi- 
tuting x p into the differential equation, 

idbe' nt = Abe mt + e iQt ej. 

Equating coefficients of the exponential, 


(idI-A)b = e j . (25) 

As id is not an eigenvalue of A, the solution may be formed as 

b = (idf - A)-%. 

For the problem of interest, the non-secular particular solution is a linear 
combination of particular solutions of this type, each scaled by the appropri- 
ate constant coefficient found in the differential equation. In the system 


u = /\u + v. 
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say that the individual contributions to the particular solution arc now writ- 
ten as 

U fli = C,b, ; e' n ''. 

where each c r is the appropriate scaling constant. Then, the full particular 
solution is the sum of these solutions, 

u p = J2 u p’- 

i 

Combining with the homogeneous solution, the complete solution for u is 

u = T ( e " n,f - eAt ) c > b >- 

i 

Note that this solution satisfies the requirement of zero initial conditions for 
the variational equations. 

A. 3. 2 Secular Solution 

Next, examine the case where the nonhomogeneous vector v contains periodic 
(or exponential) contributions with frequencies equal to one of the eigenvalues 
of A. In this case, a secular solution will result. 

Consider a general first-order nonhomogeneous linear system of differen- 
tial equations given by 

x = Ax + e Xit Bj, 

where X L is the i th eigenvector of the matrix A. 

The particular solution of such a system may again be found using the 
method of undetermined coefficients. Say that the particular solution is of 
the form 

x p = (at + b)e A C 

where the vectors a and b are constant vectors which are to be determined. 
Substituting x p into the differential equation, 

(a + A, at -fi Ajb)e 1 — /l (a/ 4- b)e A ^ -f- e A/ ^Cj. 

Equating coefficients of t and unity, 


(A i I — /l) a = 0 


( 26 ) 
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and 


(A Z I — j4)b — — a + 6j. 
Equation (26) may be solved to give 


(27) 


a = 

where x* is the i xh eigenvector of A and a is yet to be determined. Substi- 
tuting into Equation (27), 


(A J — ,4)b = -ax,. + e r 

Xext, let y 7 denote the i th left eigenvector of A, satisfying 


yT {Xi I- A) = 0 T 


(28) 


The elements of y 7 - describe the linear dependence of the rows of A J — A. 
Therefore, a summation of the weighted rows of Equation (28) allows the 
elimination of b. This is performed by writing 


0 T = y/ r (Ai/ - -4)b = -y, T cvx I + y, T ej 

= -yi T axi + t/ij. 

where y t] is the j th element of y t . Solving for cv. 

_ Vij 


a = 


y r r x, ' 


(29) 


Let the columns of \,1 — A be given by the column vectors c, . Therefore, 
an n-dimensional Equation (28) may be represented as 

[ ci c 2 ••• c„ ] b = - qx , + ej. 

Now, examine all but. the j th row of this equation: 

[ ci £2 • • • Cj, ] b = Mb 

= -ax;, 


where an under bar denotes the removal of row j. 

The eigenvalue A, is assumed to be distinct. Therefore, it is possible to 
arbitrarily sot the j lh element of b, bj, to unity, and solve for the remainder 
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of the vector. By setting bj to unity, tire contributions associated with the 
j th column of A t I — .4 may be separated from the remainder of the equation 
as 

Mb = -nx, — c n 

where an overbar denotes the removal of column j. Because tire eigenvalue is 
distinct, there is no zero eigenvalue, and a single row and column have been 
removed from A J — A , the remaining matrix, A/, is noirsingular. Therefore, 
solving for b. 

b = -A/ 1 j . (30) 

The vector b is then formed by augmenting b with the inclusion of unity as 
the 7 th element. 

To illustrate this procedure, consider the following example: 


Here, the exponential multiplier associated with the nonhomogeneous part 
of the equation is 1, which is also an eigenvalue of the system matrix. The 
corresponding eigenvector is 


xi = 


1 

1 


and the associated left eigenvector is 

yi = [ 3 -i ] T 

Equation (29) gives a as 

3 _ 3 

a_ (3)(l) + (-l)(l) “ 2’ 

Equation (30) gives the single element vector b as 


b= -0/3) 


|(1) + (-3) 


1 

2 ' 


Therefore, the particular solution is 

3 r i 

i 


x p = \ 2 


1 

1/2 



For the problem of interest, the particular solution is a linear combination 
of particular solutions of this type, each scaled by the appropriate constant 
coefficient found in the differential equation. In the system 

u - An + v, 


say that the individual contributions to the particular solution are now writ- 
ten as 

u pi — c\ (a^Xi + b 7 ) c , 

where c x is the appropriate scaling constant. Then, the full particular solution 
is the sum of these solutions, 


i 

Combining with the homogeneous solution, the complete solution for u is 

u = ~ e ' 4t 7> b - + u p- 

i 

Note that this solution satisfies the requirement of zero initial conditions for 
the variational equations. 

A. 3. 3 State Transition Matrix 

In both the secular and non-secular cases, the exponential state transition 
matrix is constructed in the usual fashion. Letting P denote the matrix of 
eigenvectors of the matrix A , and A the diagonal matrix of eigenvalues, the 
exponential matrix is given by 

e At = PAP - 1 . 


A. 4 Out-of-Plane Solution 

This approach is now applied to the z and i-components of Equation (22), 
which arc not coupled with the in-plane components. These out-of-plane 
components are 


U s - 


0 1 
-a 3 0 


IJ Z + 3A 4 


0 0 0 
2 0 X 


( 31 ) 
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where 


U z = 


For a single vector u of IJ Z , 


dz 

dz 

dz 


dx h 

dyi, 

dzh 


dz 

dz 

()z 


dx h 

dy h 

d^h - 

n 


u = 


0 1 
A z 0 


u + v, 


(32) 


where v is the appropriate nonhomogeneous column of Equation (31). 

For each non-zero column v, the periodic forcing function is of frequency 
which is the same or approximately the same as the natural frequency of the 
system. (If higher-order frequency compensation is employed, the low-order 
frequencies of the two non-zero columns will, in fact, be identical.) Therefore, 
the secular response associated with each non-zero forcing column is exam- 
ined. Additionally, the non-secular response is examined for the case where 
the third-column forcing function possesses the nearly resonant frequency of 
the linear solution. 


A. 4.1 Secular 

First consider the secular case, where v is of the form 


0 

Asin(y/A^t + 1/ ; ) 

This corresponds to the first nonhomogeneous column of Equation (31). (If 
the third column is treated as secular, it differs only by a 90° shift in phase.) 
Using exponential representation, 


v = A\ 


0 

jis/A^t 


“b At) 


0_ 

- 1 \/ A i t 


whore 


and 


A\ = — (sin v — i cos xb) 
, A, 

A 2 = — (sin u > + i cos ip ) . 
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The particular solution may now bo determined using the method of Sec- 
tion A. 3. 2, treating each of the two exponential vectors of v individually. 

Using the earlier notation, let A] = i\/A^ and A 2 = —is/A^. The eigen- 
vectors are 


*i - 

the left eigenvectors are 


1 


and X 2 - 


— iy/A^ _ 

yi = [ i'/M 1 ] T and ya = [ -i\fM 1 ] T 


determined coefficients gives 


For the nonhomogeneous contribution given by e iv ^ t e 2 , the method of un- 

1 


and 


b = 


1 


2.4,3 y/ A 3 

1 

Therefore, this part of the particular solution is 

1 i 


u p i — A\ 



i 


1 

y/Ai 

t + 

2 

1 




•2/I3 VK 
1 




In similar fashion, the part of the particular solution corresponding to 
e -i j s fd rme d. Combining gives the full particular solution as 


Up = A\ 



+ Ai 


1 

2 


V 


i 


1 i 

\ 

y/Al 

t + 

2/1 3 


1 


1 

/ 

1 

t + 

1 i 

2 .4 3 yj A3 

\ 

1 


1 

/ 
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Referring to Section A. 3. 2 , the behavior of the position partial derivative 
over 180 days is presented in Figure 10, for a pair of fundamental values of the 
phase angle. As previously mentioned, the hub position error components are 
taken to be 1 km; therefore, in this and all other figures, the contributions to 
the telescope motion variation is the same as the position partial derivatives 
themselves. 

A. 4. 2 Non-Secular 

Next, consider the response associated with Equation (32), where the vector 
v is a nearly-rosonant vector of the form 

A cos(u4 + (j))e 2 - 

From the earlier analysis, the phase angle (f> is taken as 0+9 0°. The frequency 
uj is taken as the natural linear frequency of the in-plane motion. As shown 
in the earlier analysis, this frequency is given by 

to = \J n 1 — .A ^ / 2 + y/(n 2 — A$/2) 2 — (n a + 2/la )(?v . 2 — A 3 ) . 
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For comparison, to is approximately 0.035385 rad/day, while s/A^ is approx- 
imately 0.034148 rad/day. 

Following the method of Section A. 3.1, the vector b is constructed as 

(to I — A) 1 & 2 - 

The resulting position partials in the solution vector u are presented in Fig- 
ure 11. 

A. 5 In-Plane Particular Solution 

The in-plane ( x and y ) components of Equation (22) are now treated. These 
differential equations are of the form 


0 

0 

1 

0 



0 

0 

0 ' 

0 

0 

0 

1 

u xl 

+ 3/U 

0 

0 

y 

0 

n 2 + 2A 3 

0 

0 

2 n 

-2.x 

Z 

0 

n 2 - 4, 

—2 n 

0 



y 

X 

o _ 
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where 


dx 

Ox 

dx 

dx h 

Otjh 

dT h 

dy_ 

Oy_ 

Oy_ 

dx h 

Oyh 

dzh 

Ox 

dx 

dx 

dx h 

dyi,. 

d7h 

8 y 

8 y 

Oy 

0 r h 

( hjh 

dz h 


For a single vector u of U :r , n 

0 0 10 

0 0 0 1 

U “ n 2 + 2 An 0 0 2 n 

0 n 2 - A 3 -2/7, 0 


U + V, 


(34) 


whore v is the appropriate column of the above matrix. 

From the earlier analysis of the linearized circular restricted three-bodv 
problem, the system matrix of Equation (34) has one eigenvalue correspond- 
ing to a convergent solution, one c orresponding to a divergent solution, and 
a complex conjugate pair which corresponds to a periodic solution. In that 
analysis, it was specified that the initial conditions would be selected in such 
a manner as to excite only the periodic modes. Therefore, as in the previ- 
ous section, it is assumed here that x and y in Equation (34) exhibit that 
behavior; z has already been seen to be periodic. 


A. 5.1 Secular 

First, examine the secular solutions associated with resonant forcing vectors 
v. 

Again using the notation of Section A. 3. 2, consider the eigenvalues given 
by Ai = iuj and Ao = —iuj. The corresponding eigenvectors are 



■i.2nuj 



—i2nu 

Xl = 

<N 

* 3 
1 Fi 
1 

and 

x 2 - 

-x 

—InJ 1 


-MX 



UxJX 
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where 

X — + n ~ + 2 A 3 ; 

the left eigenvectors are 


yi = 


-b 2 A3) 
(n 2 - A 3 )x 
—2 nuJ 2 
iu>X 


and 


y 2 = 


—i2nuj(n 2 + 2 . 43 ) 
(n 2 - A-i)x 
- —2nw 2 

—iojx 


Consider the contribution associated with the second column of Equa- 
tion (33) as the vector v of Equation (34). The order 1 solutions for x and 
y are of the form 


3 A 4 X = -A cos {(jjt, + d) and 3A 4 y — kA sin(u,’< + </>), 


where 


k.= 


X 

2 urn 


These functions may then again be written as combinations of complex ex- 
ponentials. 

For the nonhornogeneous contribution to v of e’ wt e^, the scalar cv is given 


by 


a = 


2 mo 2 

d 


where 


d — u/’ + 5(n~ + A%)u>^ — (5 n~ — 4A3)(r?/ 4- 2A'$)uj 2 — (n" — A- 4 ){n~ -t- 2A\p) 2 . 


The associated vector b is 


4n 2 u> 2 /d — i/u) 

2 n/( + i4nxu 3 /(dQ 
1 

— 2 nui 2 x(u 2 - ri 2 + Az)/(dQ T i2nui/£ 
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where 

(, = u) + n — A3. 

For the contribution of e” lw< e3, the scalar a is the same, and the associated 
vector b is the complex conjugate. 

For the contribution of 


and 

— 2 n/x ~ i4nu>*/d 
X 2 /d - i/u, 

2mo 2 (u> 2 — n 2 — 2A :i )/d — i‘2nuj/x 

1 

For the contribution of t lujl e^. the resulting a and b are the complex conju- 
gates of these. 

Combining into the solution vector u, the resulting position partial deriva- 
tives, in the second column of matrix U xy , are presented in Figures 12(a) 
(./-component) and 12(b) (//-component). 






Ox 

Figure 13: — — and 
dx h 



with resonance 


Next, the same procedure is followed for the first column. The solution is 
constructed using the same a scalars and b vectors as above. The resulting 
position partial derivatives are presented in Figures 13(a) (^-component) 
and 13(b) (y-component). 

Finally, the procedure is followed for the third column. If the low-order 
solution for z is taken as being resonant with the system matrix, the solution 
associated with the forcing function 3A f \ze% may be constructed using the 
same a and b expressions as above. The resulting position partial derivatives 
are presented in Figures 14(a) (^-component) and 14(b) (y-component). 

A. 5. 2 Non-Secular 

Finally, the case is treated where the uncompensated low-order solution for 2 
is considered to be non-resonant with the in-plane system. The third-column 
forcing vector is taken to be of the form 

+ <p)e 3 . 

Again following the method of Section A.3.1, the vector b is constructed as 

(uf - /l) _1 e 3 . 
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The resulting solution vector u is presented in Figures 15(a) (^-component) 
and 15(b) (y-component). 

A. 6 Summary 

From the plots of this appendix, it is recognized that all of the small- 
magnitude contributions to the solution of the variational equations are in 
the millimeter range. Over 180 days, the large- magnitude contributions reach 
hundreds of millimeters. However, even these contributions remain extremely 
small over the first 90 days from epoch. Therefore, it is concluded that the 
errors which are induced by 1 km errors in hub position are sufficiently small 
that they may be ignored. 
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B Relative Motion Near Earth-Moon L 2 Li- 
bration Point 

At the sponsor’s request, a rudimentary analysis was performed into the 
approach which would be taken in an investigation of relative motion near 
the Earth-Moon L 2 point. In this system. Earth and Moon are the pri- 
mary bodies; here, the fundamental periods in and out of the zy-plane axe 
approximately 14.650 days and 15.275 days, respectively. 

Consider the ease of a hub orbiting the /, 2 point at roughly a distance 
of 30,000 km. with a telescope located 100 m from the hub. The effects of 
various terms of the elliptical restricted expansion over 15 days are given in 
Table 5. 

Table 5: Along- Ellipsoid Effect of Earth-Moon Perturbation on Solution (15 
days) 


perturbing 

term 

effect 

(in) 

re 

21.5 

rr h 

114.1 

rr h e 

8.4 

rrh 

64.4 

rr h 3 

31.5 

?•)'/, 4 

14.9 


The effects are shown for terms which contribute to as low as approxi- 
mately 10 m over the 15 days. It is noted that the largest contributions from 
perturbations by Sun during this time period are on the order of less than 
1 rn. Perturbing forces from solar radiation pressure and spacecraft thruster 
firing are on the same order as for the Sun-Earth system. 
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C Calculation of Luminosity Reduction Due 
to Partial Eclipse 

The calculation of the luminosity reduction factor, cr, is shown here. This 
symbol is a variable within Equation (8), which is used to compute the force 
per unit mass due to solar radiation pressure. 



Figure 16: Solar Radiation Shadowing Geometry 


Figure 16 shows the shadowing geometry. In the figure, line CS is normal 
to line PQ; line CD is normal to the Sun-Eartli line, SE; line SO is normal 
to line CO. The distance Rt> is the Earth’s effective obscuration radius. Dis- 
tances Re and Rs are the mean radii of Earth and Sun, respectively. Solar 
shadowing is based on projecting 2 /?& onto line PQ. 

The calculation of a is a matter of the geometry of two intersecting circles. 
Figure 17 shows the additional geometry and the notation used in the cal- 
culation. The projection of the obscuration disk on the plane normal to the 
spacecraft-Sun line (CS) is an ellipse of negligible eccentricity. This ellipse 
is very closely approximated as a circle. The comment lines in the accom- 
panying FORTRAN code describes the notations used in Figures 16 and 17. 
The routine calculates the solar radiation force per unit mass at any distance 
greater than the Sun- Earth distance. The code also calculates the diminished 
radiation force for all spacecraft trajectories that pass through any portion 
of the Earth’s shadow. It is assumed that the radiation force falls off linearly 
in proportion to the percentage of solar area lost to obscuration. 
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PROGRAM SOLRAD 
C 

IMPLICIT REAL*8 (A-H , O-Z) 

REAL*8 LRF, MASS 
COMMON/ INPUTS/CR , AREA , MASS 
C 

C EXAMPLE DRIVER FOR SOLAR RADIATION FORCE 

C 

CR=1 . 5D0 
AREA=10 . DO 
MASS=2000 . DO 
SE= 149597870 . DO 
C SEE FIG 16: 

WRITE (* , *) > INPUT CE AND CD = ' 

READ ( * , * ) CE , CD 
C 

CALL LRFACTOR(SE,CE,CD,LRF) 

WRITE(* , *) ' LRF = ’ ,LRF 
C 

CS=DSQRT (SE**2+CE**2) 

CALL RFORCE ( LRF , CS , F) 

WRITE(* , 101) F 

101 FORMAT(' SOLAR RADIATION FORCE/MASS (M/S“2) = ’.1PD11.4) 
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c 


c 

c 

c 

c 

•c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


STOP 

END 

SUBROUTINE LRFACTOR ( SE , CE , CD , LRF ) 

IMPLICIT REAL*8(A-H,0~Z) 

REAL*8 LRF,L 
DATA RS/ 696000 . DO/ 

DATA RE/6378. 14D0/ 

DATA PI/3 . 141592653589793D0/ 

LRF = LUMINOSITY REDUCTION FACTOR (OUTPUT) 

FIG 16 CALCULATIONS: 

RS = SOLAR RADIUS (KM) 

RE = EARTH RADIUS (KM) 

SE = SUN EARTH DISTANCE (KM) (INPUT) 

CE = SATELLITE TO EARTH DISTANCE (KM) (INPUT) 

CE IS A PORTION OF CO. 

CD = NORMAL DISTANCE OF SATELLITE TO SUN-EARTH 
LINE (KM) (INPUT) (GE ZERO) 

CS = SATELLITE TO SUN DISTANCE (KM) 

CO = SATELLITE TO SUN (THRU E) OFFSET DISTANCE (KM) 

OS = PERPENDICULAR DISTANCE OF SUN FROM CO (KM) 

CO + OS + CS FORMS A RIGHT TRIANGLE WITH 
HYPOTENUSE CS . 

PS = NORTHERN LIMIT OF SOLAR OSBSCURATION ALONG LINE PQ 
NORMAL TO CS . PS MEASURED FROM SOLAR CENTER S. 

FIG 17 CALCULATIONS: 

RO = RADIUS OF OBSCURATION DISK (KM) 

C = SEPARATION BETWEEN GEOMETRIC CENTER OF SOLAR DISK 
AND GEOMETRIC CENTER OF OBSCURATION DISK (KM) 

= SMALLA + SMALLB 
= SO’ 
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C L = DISTANCE FROM LINE SEGMENT C TO EITHER POINT OF 

C -INTERSECTION OF THE 2 DISK PERIMETERS (KM) 

C L**2 = R0**2 - SMALLA**2 = Rl**2 - SMALLB**2 (RIGHT 

C TRIANGLES) 

C THETAO = DATAN (L/SMALLA) 

C THETA 1 = DATAN (L/SMALLB) 

C SO = SECTOR AREA SUBTENDED BY 2*THETA0 (KM**2) 

C SI = SECTOR AREA SUBTENDED BY 2*THETA1 (KM**2) 

C ASHADOW = AREA OF INTERSECTION OF OBSCURATION DISK AND 

C SOLAR DISK (KM**2) 

C SUNAREA = AREA OF SOLAR DISK (KM**2) 

C 

C SEE FIG 16 FOR THE FOLLOWING: 

C 

ED=DSQRT (CE**2-CD**2) 

C COSINE(LAMBDA) 

COSL=ED/CE 
CO=SE*COSL+CE 
C SINE(LAMBDA) 

SINL=CD/CE 

OS=SE*SINL 

CA=DSQRT (CE**2-RE**2) 

CS=DSQRT ( (SE+ED) **2+CD**2) 

C OBSCURATION OCCURS WHEN DABS (PS) .LE.RS 

PS=CS* (RE*CO-OS*CA) / (CA*CO+OS*RE) 

IF(PS.LT.-RS) THEN 
C NO ECLIPSE 

LRF=1 . ODO 
RETURN 
ENDIF 

IF(PS.GT.RS) THEN 
C TOTAL ECLIPSE 

LRF=0 . ODO 
RETURN 
ENDIF 
C 

C PARTIAL ECLIPSE (DABS(PS) .LE.RS) 

C 
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C COSINE (BETA - ALPHA) 

COSBA=CO/CS 

C "NORTHERN" RADIUS OF OBSCURATION MEASURED FROM O’ 

C "NORTH" TO P ALONG LINE PQ: 

SOP=OS/COSBA 

RADIUSN=PS+SOP 

RO=RADIUSN 

C 

C SEE FIG 17 FOR THE FOLLOWING: 

C 

C=SOP 

SMALLB= (RS**2-R0**2+C**2) / (2 . DOC) 

SMALLA=C-SMALLB 
L=DSQRT (RS**2-SMALLB**2) 

THETAO=DATAN (L/SMALLA) 

THETA1=DATAN (L/SMALLB) 

SO=RO**2*THETAO 
S1=RS**2*THETA1 
ASHADOW=SO+Sl-C*L 
SUNAREA=PI*RS**2 
LRF=1 . DO-ASHADOW/SUNAREA 
RETURN 
END 
C 

SUBROUTINE RFORCE (LRF , CS , F) 

IMPLICIT REAL+8 (A-H , O-Z) 

REAL *8 LRF .MASS 
COMMON/INPUTS/CR , AREA , MASS 
C 

C LRF = LUMINOSITY REDUCTION FACTOR (INPUT) 

C CS = HELIOCENTRIC DISTANCE TO NGST (KM) (INPUT) 

C F = SOLAR RADIATION FORCE PER UNIT MASS 

C (NEWTONS/KG = M/S~2) (OUTPUT) 

C 

C CR = SOLAR FLUX REFLECTION PARAMETER 0 <= CR <= 2 

C (INPUT THRU COMMON) 

C AREA = NGST AREA PROJECTED NORMAL TO SUN LINE CE 

C (METERS**2) 
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D Constant Parameters 


Table 6: Physical Constants [5] 


gravitational parameter, 
Sun alone 

M i 

132,712,440,017.987 km 3 /s 3 

gravitational parameter, 
Earth alone 

L L 2 

398,600.4415 km 3 /s 2 

gravitational parameter, 
Earth-Moon 

A* 2 

403,503.236 km 3 /s 2 

gravitational parameter, 
Moon alone 

Ms 

4,902.8003 km 3 /s 2 

astronomical unit 

AU 

149,597,870.691 km 

mean Earth-Moon baryccn- 
ter distance from Sun 


1.000001018 AU 

eccentricity of Earth-Moon 
baryccnter orbit about Sun 

e 

0.01670862 

mean motion of Earth-Moon 
barycenter orbit about Sun 

n 

0.199106385 xl0“ 6 rad/s 

L ‘2 distance ratio 

7 

0.01007824 


Table 7: Computed Values and Coefficients 


reference Sun- 
Earth distance 

D 

149,618,905.218739 km 

gravitational 

A, 

1.16556055765939 x 10“ 3 1/day 2 

coefficients 

A. i 

5.84525170441422 x lO -10 l/(km-day 2 ) 

(see page 126) 

A & 

3.86396147244215 x 10- le l/(km 2 -day 2 ) 



2.56240418152728 x 10“ 22 l/(km 3 -day 2 ) 
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